Actual source code: normm.c


  2: #include <petsc/private/matimpl.h>

  4: typedef struct {
  5:   Mat         A;
  6:   Vec         w,left,right,leftwork,rightwork;
  7:   PetscScalar scale;
  8: } Mat_Normal;

 10: PetscErrorCode MatScale_Normal(Mat inA,PetscScalar scale)
 11: {
 12:   Mat_Normal *a = (Mat_Normal*)inA->data;

 15:   a->scale *= scale;
 16:   return(0);
 17: }

 19: PetscErrorCode MatDiagonalScale_Normal(Mat inA,Vec left,Vec right)
 20: {
 21:   Mat_Normal     *a = (Mat_Normal*)inA->data;

 25:   if (left) {
 26:     if (!a->left) {
 27:       VecDuplicate(left,&a->left);
 28:       VecCopy(left,a->left);
 29:     } else {
 30:       VecPointwiseMult(a->left,left,a->left);
 31:     }
 32:   }
 33:   if (right) {
 34:     if (!a->right) {
 35:       VecDuplicate(right,&a->right);
 36:       VecCopy(right,a->right);
 37:     } else {
 38:       VecPointwiseMult(a->right,right,a->right);
 39:     }
 40:   }
 41:   return(0);
 42: }

 44: PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
 45: {
 46:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 48:   Vec            in;

 51:   in = x;
 52:   if (Na->right) {
 53:     if (!Na->rightwork) {
 54:       VecDuplicate(Na->right,&Na->rightwork);
 55:     }
 56:     VecPointwiseMult(Na->rightwork,Na->right,in);
 57:     in   = Na->rightwork;
 58:   }
 59:   MatMult(Na->A,in,Na->w);
 60:   MatMultTranspose(Na->A,Na->w,y);
 61:   if (Na->left) {
 62:     VecPointwiseMult(y,Na->left,y);
 63:   }
 64:   VecScale(y,Na->scale);
 65:   return(0);
 66: }

 68: PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
 69: {
 70:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 72:   Vec            in;

 75:   in = v1;
 76:   if (Na->right) {
 77:     if (!Na->rightwork) {
 78:       VecDuplicate(Na->right,&Na->rightwork);
 79:     }
 80:     VecPointwiseMult(Na->rightwork,Na->right,in);
 81:     in   = Na->rightwork;
 82:   }
 83:   MatMult(Na->A,in,Na->w);
 84:   VecScale(Na->w,Na->scale);
 85:   if (Na->left) {
 86:     MatMultTranspose(Na->A,Na->w,v3);
 87:     VecPointwiseMult(v3,Na->left,v3);
 88:     VecAXPY(v3,1.0,v2);
 89:   } else {
 90:     MatMultTransposeAdd(Na->A,Na->w,v2,v3);
 91:   }
 92:   return(0);
 93: }

 95: PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y)
 96: {
 97:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 99:   Vec            in;

102:   in = x;
103:   if (Na->left) {
104:     if (!Na->leftwork) {
105:       VecDuplicate(Na->left,&Na->leftwork);
106:     }
107:     VecPointwiseMult(Na->leftwork,Na->left,in);
108:     in   = Na->leftwork;
109:   }
110:   MatMult(Na->A,in,Na->w);
111:   MatMultTranspose(Na->A,Na->w,y);
112:   if (Na->right) {
113:     VecPointwiseMult(y,Na->right,y);
114:   }
115:   VecScale(y,Na->scale);
116:   return(0);
117: }

119: PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
120: {
121:   Mat_Normal     *Na = (Mat_Normal*)N->data;
123:   Vec            in;

126:   in = v1;
127:   if (Na->left) {
128:     if (!Na->leftwork) {
129:       VecDuplicate(Na->left,&Na->leftwork);
130:     }
131:     VecPointwiseMult(Na->leftwork,Na->left,in);
132:     in   = Na->leftwork;
133:   }
134:   MatMult(Na->A,in,Na->w);
135:   VecScale(Na->w,Na->scale);
136:   if (Na->right) {
137:     MatMultTranspose(Na->A,Na->w,v3);
138:     VecPointwiseMult(v3,Na->right,v3);
139:     VecAXPY(v3,1.0,v2);
140:   } else {
141:     MatMultTransposeAdd(Na->A,Na->w,v2,v3);
142:   }
143:   return(0);
144: }

146: PetscErrorCode MatDestroy_Normal(Mat N)
147: {
148:   Mat_Normal     *Na = (Mat_Normal*)N->data;

152:   MatDestroy(&Na->A);
153:   VecDestroy(&Na->w);
154:   VecDestroy(&Na->left);
155:   VecDestroy(&Na->right);
156:   VecDestroy(&Na->leftwork);
157:   VecDestroy(&Na->rightwork);
158:   PetscFree(N->data);
159:   return(0);
160: }

162: /*
163:       Slow, nonscalable version
164: */
165: PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
166: {
167:   Mat_Normal        *Na = (Mat_Normal*)N->data;
168:   Mat               A   = Na->A;
169:   PetscErrorCode    ierr;
170:   PetscInt          i,j,rstart,rend,nnz;
171:   const PetscInt    *cols;
172:   PetscScalar       *diag,*work,*values;
173:   const PetscScalar *mvalues;

176:   PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
177:   PetscArrayzero(work,A->cmap->N);
178:   MatGetOwnershipRange(A,&rstart,&rend);
179:   for (i=rstart; i<rend; i++) {
180:     MatGetRow(A,i,&nnz,&cols,&mvalues);
181:     for (j=0; j<nnz; j++) {
182:       work[cols[j]] += mvalues[j]*mvalues[j];
183:     }
184:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
185:   }
186:   MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
187:   rstart = N->cmap->rstart;
188:   rend   = N->cmap->rend;
189:   VecGetArray(v,&values);
190:   PetscArraycpy(values,diag+rstart,rend-rstart);
191:   VecRestoreArray(v,&values);
192:   PetscFree2(diag,work);
193:   VecScale(v,Na->scale);
194:   return(0);
195: }

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

200:    Collective on Mat

202:    Input Parameter:
203: .   A  - the (possibly rectangular) matrix

205:    Output Parameter:
206: .   N - the matrix that represents A'*A

208:    Level: intermediate

210:    Notes:
211:     The product A'*A is NOT actually formed! Rather the new matrix
212:           object performs the matrix-vector product by first multiplying by
213:           A and then A'
214: @*/
215: PetscErrorCode  MatCreateNormal(Mat A,Mat *N)
216: {
218:   PetscInt       n,nn;
219:   Mat_Normal     *Na;
220:   VecType        vtype;

223:   MatGetSize(A,NULL,&nn);
224:   MatGetLocalSize(A,NULL,&n);
225:   MatCreate(PetscObjectComm((PetscObject)A),N);
226:   MatSetSizes(*N,n,n,nn,nn);
227:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);

229:   PetscNewLog(*N,&Na);
230:   (*N)->data = (void*) Na;
231:   PetscObjectReference((PetscObject)A);
232:   Na->A      = A;
233:   Na->scale  = 1.0;

235:   MatCreateVecs(A,NULL,&Na->w);

237:   (*N)->ops->destroy          = MatDestroy_Normal;
238:   (*N)->ops->mult             = MatMult_Normal;
239:   (*N)->ops->multtranspose    = MatMultTranspose_Normal;
240:   (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal;
241:   (*N)->ops->multadd          = MatMultAdd_Normal;
242:   (*N)->ops->getdiagonal      = MatGetDiagonal_Normal;
243:   (*N)->ops->scale            = MatScale_Normal;
244:   (*N)->ops->diagonalscale    = MatDiagonalScale_Normal;
245:   (*N)->assembled             = PETSC_TRUE;
246:   (*N)->preallocated          = PETSC_TRUE;

248:   MatSetOption(*N,MAT_SYMMETRIC,PETSC_TRUE);
249:   MatGetVecType(A,&vtype);
250:   MatSetVecType(*N,vtype);
251: #if defined(PETSC_HAVE_DEVICE)
252:   MatBindToCPU(*N,A->boundtocpu);
253: #endif
254:   return(0);
255: }