Actual source code: normm.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/matimpl.h

      5: typedef struct {
      6:   Mat A;
      7:   Vec w;
      8: } Mat_Normal;

     12: PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
     13: {
     14:   Mat_Normal     *Na = (Mat_Normal*)N->data;

     18:   MatMult(Na->A,x,Na->w);
     19:   MatMultTranspose(Na->A,Na->w,y);
     20:   return(0);
     21: }

     25: PetscErrorCode MatDestroy_Normal(Mat N)
     26: {
     27:   Mat_Normal     *Na = (Mat_Normal*)N->data;

     31:   PetscObjectDereference((PetscObject)Na->A);
     32:   VecDestroy(Na->w);
     33:   PetscFree(Na);
     34:   return(0);
     35: }
     36:
     37: /*
     38:       Slow, nonscalable version
     39: */
     42: PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
     43: {
     44:   Mat_Normal        *Na = (Mat_Normal*)N->data;
     45:   Mat               A = Na->A;
     46:   PetscErrorCode    ierr;
     47:   PetscInt          i,j,rstart,rend,nnz;
     48:   const PetscInt    *cols;
     49:   PetscScalar       *diag,*work,*values;
     50:   const PetscScalar *mvalues;
     51:   PetscMap          cmap;

     54:   PetscMalloc(2*A->N*sizeof(PetscScalar),&diag);
     55:   work = diag + A->N;
     56:   PetscMemzero(work,A->N*sizeof(PetscScalar));
     57:   MatGetOwnershipRange(A,&rstart,&rend);
     58:   for (i=rstart; i<rend; i++) {
     59:     MatGetRow(A,i,&nnz,&cols,&mvalues);
     60:     for (j=0; j<nnz; j++) {
     61:       work[cols[j]] += mvalues[j]*mvalues[j];
     62:     }
     63:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
     64:   }
     65:   MPI_Allreduce(work,diag,A->N,MPIU_SCALAR,MPI_SUM,N->comm);
     66:   MatGetPetscMaps(A,PETSC_NULL,&cmap);
     67:   PetscMapGetLocalRange(cmap,&rstart,&rend);
     68:   VecGetArray(v,&values);
     69:   PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));
     70:   VecRestoreArray(v,&values);
     71:   PetscFree(diag);
     72:   return(0);
     73: }

     77: /*@
     78:       MatCreateNormal - Creates a new matrix object that behaves like A'*A.

     80:    Collective on Mat

     82:    Input Parameter:
     83: .   A  - the (possibly rectangular) matrix

     85:    Output Parameter:
     86: .   N - the matrix that represents A'*A

     88:    Level: intermediate

     90:    Notes: The product A'*A is NOT actually formed! Rather the new matrix
     91:           object performs the matrix-vector product by first multiplying by
     92:           A and then A'
     93: @*/
     94: PetscErrorCode PETSCMAT_DLLEXPORT MatCreateNormal(Mat A,Mat *N)
     95: {
     97:   PetscInt       m,n;
     98:   Mat_Normal     *Na;

    101:   MatGetLocalSize(A,&m,&n);
    102:   MatCreate(A->comm,N);
    103:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
    104:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);
    105:
    106:   PetscNew(Mat_Normal,&Na);
    107:   Na->A     = A;
    108:   PetscObjectReference((PetscObject)A);
    109:   (*N)->data = (void*) Na;

    111:   VecCreateMPI(A->comm,m,PETSC_DECIDE,&Na->w);
    112:   (*N)->ops->destroy     = MatDestroy_Normal;
    113:   (*N)->ops->mult        = MatMult_Normal;
    114:   (*N)->ops->getdiagonal = MatGetDiagonal_Normal;
    115:   (*N)->assembled        = PETSC_TRUE;
    116:   (*N)->N                = A->N;
    117:   (*N)->M                = A->N;
    118:   (*N)->n                = A->n;
    119:   (*N)->m                = A->n;
    120:   return(0);
    121: }