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 MatIncreaseOverlap_Normal(Mat A,PetscInt is_max,IS is[],PetscInt ov)
 45: {
 46:   Mat_Normal     *a = (Mat_Normal*)A->data;
 47:   Mat            pattern;

 51:   if (ov < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap specified");
 52:   MatProductCreate(a->A,a->A,NULL,&pattern);
 53:   MatProductSetType(pattern,MATPRODUCT_AtB);
 54:   MatProductSetFromOptions(pattern);
 55:   MatProductSymbolic(pattern);
 56:   MatIncreaseOverlap(pattern,is_max,is,ov);
 57:   MatDestroy(&pattern);
 58:   return(0);
 59: }

 61: PetscErrorCode MatCreateSubMatrices_Normal(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[])
 62: {
 63:   Mat_Normal     *a = (Mat_Normal*)mat->data;
 64:   Mat            B = a->A, *suba;
 65:   IS             *row;
 66:   PetscInt       M;

 70:   if (a->left || a->right || irow != icol) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Not implemented");
 71:   if (scall != MAT_REUSE_MATRIX) {
 72:     PetscCalloc1(n,submat);
 73:   }
 74:   MatGetSize(B,&M,NULL);
 75:   PetscMalloc1(n,&row);
 76:   ISCreateStride(PETSC_COMM_SELF,M,0,1,&row[0]);
 77:   ISSetIdentity(row[0]);
 78:   for (M = 1; M < n; ++M) row[M] = row[0];
 79:   MatCreateSubMatrices(B,n,row,icol,MAT_INITIAL_MATRIX,&suba);
 80:   for (M = 0; M < n; ++M) {
 81:     MatCreateNormal(suba[M],*submat+M);
 82:     ((Mat_Normal*)(*submat)[M]->data)->scale = a->scale;
 83:   }
 84:   ISDestroy(&row[0]);
 85:   PetscFree(row);
 86:   MatDestroySubMatrices(n,&suba);
 87:   return(0);
 88: }

 90: PetscErrorCode MatPermute_Normal(Mat A,IS rowp,IS colp,Mat *B)
 91: {
 92:   Mat_Normal     *a = (Mat_Normal*)A->data;
 93:   Mat            C,Aa = a->A;
 94:   IS             row;

 98:   if (rowp != colp) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Row permutation and column permutation must be the same");
 99:   ISCreateStride(PetscObjectComm((PetscObject)Aa),Aa->rmap->n,Aa->rmap->rstart,1,&row);
100:   ISSetIdentity(row);
101:   MatPermute(Aa,row,colp,&C);
102:   ISDestroy(&row);
103:   MatCreateNormal(C,B);
104:   MatDestroy(&C);
105:   return(0);
106: }

108: PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
109: {
110:   Mat_Normal     *a = (Mat_Normal*)A->data;
111:   Mat            C;

115:   if (a->left || a->right) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
116:   MatDuplicate(a->A,op,&C);
117:   MatCreateNormal(C,B);
118:   MatDestroy(&C);
119:   if (op == MAT_COPY_VALUES) ((Mat_Normal*)(*B)->data)->scale = a->scale;
120:   return(0);
121: }

123: PetscErrorCode MatCopy_Normal(Mat A,Mat B,MatStructure str)
124: {
125:   Mat_Normal     *a = (Mat_Normal*)A->data,*b = (Mat_Normal*)B->data;

129:   if (a->left || a->right) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented");
130:   MatCopy(a->A,b->A,str);
131:   b->scale = a->scale;
132:   VecDestroy(&b->left);
133:   VecDestroy(&b->right);
134:   VecDestroy(&b->leftwork);
135:   VecDestroy(&b->rightwork);
136:   return(0);
137: }

139: PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
140: {
141:   Mat_Normal     *Na = (Mat_Normal*)N->data;
143:   Vec            in;

146:   in = x;
147:   if (Na->right) {
148:     if (!Na->rightwork) {
149:       VecDuplicate(Na->right,&Na->rightwork);
150:     }
151:     VecPointwiseMult(Na->rightwork,Na->right,in);
152:     in   = Na->rightwork;
153:   }
154:   MatMult(Na->A,in,Na->w);
155:   MatMultTranspose(Na->A,Na->w,y);
156:   if (Na->left) {
157:     VecPointwiseMult(y,Na->left,y);
158:   }
159:   VecScale(y,Na->scale);
160:   return(0);
161: }

163: PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
164: {
165:   Mat_Normal     *Na = (Mat_Normal*)N->data;
167:   Vec            in;

170:   in = v1;
171:   if (Na->right) {
172:     if (!Na->rightwork) {
173:       VecDuplicate(Na->right,&Na->rightwork);
174:     }
175:     VecPointwiseMult(Na->rightwork,Na->right,in);
176:     in   = Na->rightwork;
177:   }
178:   MatMult(Na->A,in,Na->w);
179:   VecScale(Na->w,Na->scale);
180:   if (Na->left) {
181:     MatMultTranspose(Na->A,Na->w,v3);
182:     VecPointwiseMult(v3,Na->left,v3);
183:     VecAXPY(v3,1.0,v2);
184:   } else {
185:     MatMultTransposeAdd(Na->A,Na->w,v2,v3);
186:   }
187:   return(0);
188: }

190: PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y)
191: {
192:   Mat_Normal     *Na = (Mat_Normal*)N->data;
194:   Vec            in;

197:   in = x;
198:   if (Na->left) {
199:     if (!Na->leftwork) {
200:       VecDuplicate(Na->left,&Na->leftwork);
201:     }
202:     VecPointwiseMult(Na->leftwork,Na->left,in);
203:     in   = Na->leftwork;
204:   }
205:   MatMult(Na->A,in,Na->w);
206:   MatMultTranspose(Na->A,Na->w,y);
207:   if (Na->right) {
208:     VecPointwiseMult(y,Na->right,y);
209:   }
210:   VecScale(y,Na->scale);
211:   return(0);
212: }

214: PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
215: {
216:   Mat_Normal     *Na = (Mat_Normal*)N->data;
218:   Vec            in;

221:   in = v1;
222:   if (Na->left) {
223:     if (!Na->leftwork) {
224:       VecDuplicate(Na->left,&Na->leftwork);
225:     }
226:     VecPointwiseMult(Na->leftwork,Na->left,in);
227:     in   = Na->leftwork;
228:   }
229:   MatMult(Na->A,in,Na->w);
230:   VecScale(Na->w,Na->scale);
231:   if (Na->right) {
232:     MatMultTranspose(Na->A,Na->w,v3);
233:     VecPointwiseMult(v3,Na->right,v3);
234:     VecAXPY(v3,1.0,v2);
235:   } else {
236:     MatMultTransposeAdd(Na->A,Na->w,v2,v3);
237:   }
238:   return(0);
239: }

241: PetscErrorCode MatDestroy_Normal(Mat N)
242: {
243:   Mat_Normal     *Na = (Mat_Normal*)N->data;

247:   MatDestroy(&Na->A);
248:   VecDestroy(&Na->w);
249:   VecDestroy(&Na->left);
250:   VecDestroy(&Na->right);
251:   VecDestroy(&Na->leftwork);
252:   VecDestroy(&Na->rightwork);
253:   PetscFree(N->data);
254:   PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMat_C",NULL);
255:   PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_seqaij_C",NULL);
256:   PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_mpiaij_C",NULL);
257:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_seqdense_C",NULL);
258:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_mpidense_C",NULL);
259:   PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_dense_C",NULL);
260:   return(0);
261: }

263: /*
264:       Slow, nonscalable version
265: */
266: PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
267: {
268:   Mat_Normal        *Na = (Mat_Normal*)N->data;
269:   Mat               A   = Na->A;
270:   PetscErrorCode    ierr;
271:   PetscInt          i,j,rstart,rend,nnz;
272:   const PetscInt    *cols;
273:   PetscScalar       *diag,*work,*values;
274:   const PetscScalar *mvalues;

277:   PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
278:   PetscArrayzero(work,A->cmap->N);
279:   MatGetOwnershipRange(A,&rstart,&rend);
280:   for (i=rstart; i<rend; i++) {
281:     MatGetRow(A,i,&nnz,&cols,&mvalues);
282:     for (j=0; j<nnz; j++) {
283:       work[cols[j]] += mvalues[j]*mvalues[j];
284:     }
285:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
286:   }
287:   MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
288:   rstart = N->cmap->rstart;
289:   rend   = N->cmap->rend;
290:   VecGetArray(v,&values);
291:   PetscArraycpy(values,diag+rstart,rend-rstart);
292:   VecRestoreArray(v,&values);
293:   PetscFree2(diag,work);
294:   VecScale(v,Na->scale);
295:   return(0);
296: }

298: PetscErrorCode MatNormalGetMat_Normal(Mat A,Mat *M)
299: {
300:   Mat_Normal *Aa = (Mat_Normal*)A->data;

303:   *M = Aa->A;
304:   return(0);
305: }

307: /*@
308:       MatNormalGetMat - Gets the Mat object stored inside a MATNORMAL

310:    Logically collective on Mat

312:    Input Parameter:
313: .   A  - the MATNORMAL matrix

315:    Output Parameter:
316: .   M - the matrix object stored inside A

318:    Level: intermediate

320: .seealso: MatCreateNormal()

322: @*/
323: PetscErrorCode MatNormalGetMat(Mat A,Mat *M)
324: {

331:   PetscUseMethod(A,"MatNormalGetMat_C",(Mat,Mat*),(A,M));
332:   return(0);
333: }

335: PetscErrorCode MatConvert_Normal_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
336: {
337:   Mat_Normal     *Aa = (Mat_Normal*)A->data;
338:   Mat            B;
339:   PetscInt       m,n,M,N;

343:   MatGetSize(A,&M,&N);
344:   MatGetLocalSize(A,&m,&n);
345:   if (reuse == MAT_REUSE_MATRIX) {
346:     B = *newmat;
347:     MatProductReplaceMats(Aa->A,Aa->A,NULL,B);
348:   } else {
349:     MatProductCreate(Aa->A,Aa->A,NULL,&B);
350:     MatProductSetType(B,MATPRODUCT_AtB);
351:     MatProductSetFromOptions(B);
352:     MatProductSymbolic(B);
353:     MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
354:   }
355:   MatProductNumeric(B);
356:   if (reuse == MAT_INPLACE_MATRIX) {
357:     MatHeaderReplace(A,&B);
358:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
359:   MatConvert(*newmat,MATAIJ,MAT_INPLACE_MATRIX,newmat);
360:   return(0);
361: }

363: typedef struct {
364:   Mat work[2];
365: } Normal_Dense;

367: PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
368: {
369:   Mat            A,B;
370:   Normal_Dense   *contents;
371:   Mat_Normal     *a;
372:   PetscScalar    *array;

376:   MatCheckProduct(C,3);
377:   A = C->product->A;
378:   a = (Mat_Normal*)A->data;
379:   B = C->product->B;
380:   contents = (Normal_Dense*)C->product->data;
381:   if (!contents) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
382:   if (a->right) {
383:     MatCopy(B,C,SAME_NONZERO_PATTERN);
384:     MatDiagonalScale(C,a->right,NULL);
385:   }
386:   MatProductNumeric(contents->work[0]);
387:   MatDenseGetArrayWrite(C,&array);
388:   MatDensePlaceArray(contents->work[1],array);
389:   MatProductNumeric(contents->work[1]);
390:   MatDenseRestoreArrayWrite(C,&array);
391:   MatDenseResetArray(contents->work[1]);
392:   MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
393:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
394:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
395:   MatScale(C,a->scale);
396:   return(0);
397: }

399: PetscErrorCode MatNormal_DenseDestroy(void *ctx)
400: {
401:   Normal_Dense   *contents = (Normal_Dense*)ctx;

405:   MatDestroy(contents->work);
406:   MatDestroy(contents->work+1);
407:   PetscFree(contents);
408:   return(0);
409: }

411: PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
412: {
413:   Mat            A,B;
414:   Normal_Dense   *contents = NULL;
415:   Mat_Normal     *a;
416:   PetscScalar    *array;
417:   PetscInt       n,N,m,M;

421:   MatCheckProduct(C,4);
422:   if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
423:   A = C->product->A;
424:   a = (Mat_Normal*)A->data;
425:   if (a->left) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Not implemented");
426:   B = C->product->B;
427:   MatGetLocalSize(C,&m,&n);
428:   MatGetSize(C,&M,&N);
429:   if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
430:     MatGetLocalSize(B,NULL,&n);
431:     MatGetSize(B,NULL,&N);
432:     MatGetLocalSize(A,&m,NULL);
433:     MatGetSize(A,&M,NULL);
434:     MatSetSizes(C,m,n,M,N);
435:   }
436:   MatSetType(C,((PetscObject)B)->type_name);
437:   MatSetUp(C);
438:   PetscNew(&contents);
439:   C->product->data = contents;
440:   C->product->destroy = MatNormal_DenseDestroy;
441:   if (a->right) {
442:     MatProductCreate(a->A,C,NULL,contents->work);
443:   } else {
444:     MatProductCreate(a->A,B,NULL,contents->work);
445:   }
446:   MatProductSetType(contents->work[0],MATPRODUCT_AB);
447:   MatProductSetFromOptions(contents->work[0]);
448:   MatProductSymbolic(contents->work[0]);
449:   MatProductCreate(a->A,contents->work[0],NULL,contents->work+1);
450:   MatProductSetType(contents->work[1],MATPRODUCT_AtB);
451:   MatProductSetFromOptions(contents->work[1]);
452:   MatProductSymbolic(contents->work[1]);
453:   MatDenseGetArrayWrite(C,&array);
454:   MatSeqDenseSetPreallocation(contents->work[1],array);
455:   MatMPIDenseSetPreallocation(contents->work[1],array);
456:   MatDenseRestoreArrayWrite(C,&array);
457:   C->ops->productnumeric = MatProductNumeric_Normal_Dense;
458:   return(0);
459: }

461: PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
462: {
464:   C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
465:   return(0);
466: }

468: PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
469: {
470:   Mat_Product    *product = C->product;

474:   if (product->type == MATPRODUCT_AB) {
475:     MatProductSetFromOptions_Normal_Dense_AB(C);
476:   }
477:   return(0);
478: }

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

483:    Collective on Mat

485:    Input Parameter:
486: .   A  - the (possibly rectangular) matrix

488:    Output Parameter:
489: .   N - the matrix that represents A'*A

491:    Level: intermediate

493:    Notes:
494:     The product A'*A is NOT actually formed! Rather the new matrix
495:           object performs the matrix-vector product by first multiplying by
496:           A and then A'
497: @*/
498: PetscErrorCode  MatCreateNormal(Mat A,Mat *N)
499: {
501:   PetscInt       n,nn;
502:   Mat_Normal     *Na;
503:   VecType        vtype;

506:   MatGetSize(A,NULL,&nn);
507:   MatGetLocalSize(A,NULL,&n);
508:   MatCreate(PetscObjectComm((PetscObject)A),N);
509:   MatSetSizes(*N,n,n,nn,nn);
510:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);
511:   PetscLayoutReference(A->cmap,&(*N)->rmap);
512:   PetscLayoutReference(A->cmap,&(*N)->cmap);

514:   PetscNewLog(*N,&Na);
515:   (*N)->data = (void*) Na;
516:   PetscObjectReference((PetscObject)A);
517:   Na->A      = A;
518:   Na->scale  = 1.0;

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

522:   (*N)->ops->destroy           = MatDestroy_Normal;
523:   (*N)->ops->mult              = MatMult_Normal;
524:   (*N)->ops->multtranspose     = MatMultTranspose_Normal;
525:   (*N)->ops->multtransposeadd  = MatMultTransposeAdd_Normal;
526:   (*N)->ops->multadd           = MatMultAdd_Normal;
527:   (*N)->ops->getdiagonal       = MatGetDiagonal_Normal;
528:   (*N)->ops->scale             = MatScale_Normal;
529:   (*N)->ops->diagonalscale     = MatDiagonalScale_Normal;
530:   (*N)->ops->increaseoverlap   = MatIncreaseOverlap_Normal;
531:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
532:   (*N)->ops->permute           = MatPermute_Normal;
533:   (*N)->ops->duplicate         = MatDuplicate_Normal;
534:   (*N)->ops->copy              = MatCopy_Normal;
535:   (*N)->assembled              = PETSC_TRUE;
536:   (*N)->preallocated           = PETSC_TRUE;

538:   PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMat_C",MatNormalGetMat_Normal);
539:   PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_seqaij_C",MatConvert_Normal_AIJ);
540:   PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_mpiaij_C",MatConvert_Normal_AIJ);
541:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_seqdense_C",MatProductSetFromOptions_Normal_Dense);
542:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_mpidense_C",MatProductSetFromOptions_Normal_Dense);
543:   PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_dense_C",MatProductSetFromOptions_Normal_Dense);
544:   MatSetOption(*N,MAT_SYMMETRIC,PETSC_TRUE);
545:   MatGetVecType(A,&vtype);
546:   MatSetVecType(*N,vtype);
547: #if defined(PETSC_HAVE_DEVICE)
548:   MatBindToCPU(*N,A->boundtocpu);
549: #endif
550:   return(0);
551: }