1: #define PETSCMAT_DLL

 3:  #include include/private/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: }
 22: 
 25: PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
 26: {
 27:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 29: 
 31:   MatMult(Na->A,v1,Na->w);
 32:   MatMultTransposeAdd(Na->A,Na->w,v2,v3);
 33:   return(0);
 34: }

 38: PetscErrorCode MatDestroy_Normal(Mat N)
 39: {
 40:   Mat_Normal     *Na = (Mat_Normal*)N->data;

 44:   PetscObjectDereference((PetscObject)Na->A);
 45:   VecDestroy(Na->w);
 46:   PetscFree(Na);
 47:   return(0);
 48: }
 49: 
 50: /*
 51:       Slow, nonscalable version
 52: */
 55: PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
 56: {
 57:   Mat_Normal        *Na = (Mat_Normal*)N->data;
 58:   Mat               A = Na->A;
 59:   PetscErrorCode    ierr;
 60:   PetscInt          i,j,rstart,rend,nnz;
 61:   const PetscInt    *cols;
 62:   PetscScalar       *diag,*work,*values;
 63:   const PetscScalar *mvalues;

 66:   PetscMalloc(2*A->cmap.N*sizeof(PetscScalar),&diag);
 67:   work = diag + A->cmap.N;
 68:   PetscMemzero(work,A->cmap.N*sizeof(PetscScalar));
 69:   MatGetOwnershipRange(A,&rstart,&rend);
 70:   for (i=rstart; i<rend; i++) {
 71:     MatGetRow(A,i,&nnz,&cols,&mvalues);
 72:     for (j=0; j<nnz; j++) {
 73:       work[cols[j]] += mvalues[j]*mvalues[j];
 74:     }
 75:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
 76:   }
 77:   MPI_Allreduce(work,diag,A->cmap.N,MPIU_SCALAR,MPI_SUM,N->comm);
 78:   rstart = N->cmap.rstart;
 79:   rend   = N->cmap.rend;
 80:   VecGetArray(v,&values);
 81:   PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));
 82:   VecRestoreArray(v,&values);
 83:   PetscFree(diag);
 84:   return(0);
 85: }

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

 92:    Collective on Mat

 94:    Input Parameter:
 95: .   A  - the (possibly rectangular) matrix

 97:    Output Parameter:
 98: .   N - the matrix that represents A'*A

100:    Level: intermediate

102:    Notes: The product A'*A is NOT actually formed! Rather the new matrix
103:           object performs the matrix-vector product by first multiplying by
104:           A and then A'
105: @*/
106: PetscErrorCode  MatCreateNormal(Mat A,Mat *N)
107: {
109:   PetscInt       m,n;
110:   Mat_Normal     *Na;

113:   MatGetLocalSize(A,&m,&n);
114:   MatCreate(A->comm,N);
115:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
116:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);
117: 
118:   PetscNew(Mat_Normal,&Na);
119:   Na->A     = A;
120:   PetscObjectReference((PetscObject)A);
121:   (*N)->data = (void*) Na;

123:   VecCreateMPI(A->comm,m,PETSC_DECIDE,&Na->w);
124:   (*N)->ops->destroy     = MatDestroy_Normal;
125:   (*N)->ops->mult        = MatMult_Normal;
126:   (*N)->ops->multadd     = MatMultAdd_Normal;
127:   (*N)->ops->getdiagonal = MatGetDiagonal_Normal;
128:   (*N)->assembled        = PETSC_TRUE;
129:   (*N)->cmap.N           = A->cmap.N;
130:   (*N)->rmap.N           = A->cmap.N;
131:   (*N)->cmap.n           = A->cmap.n;
132:   (*N)->rmap.n           = A->cmap.n;
133:   return(0);
134: }