Actual source code: normalmat.c

  1: #include "normalmat.h"                /*I  "mat.h"  I*/

  5: /*@C
  6:    MatCreateADA - Creates a matrix M=A^T D1 A + D2.

  8:    Collective on matrix

 10:    Input Parameters:
 11: +  mat - matrix of arbitrary type
 12: .  D1 - A vector
 13: -  D2 - A vector

 15:    Output Parameters:
 16: .  J - New matrix whose operations are defined in terms of mat, D1, and D2.

 18:    Notes: 
 19:    The user provides the input data and is responsible for destroying
 20:    this data after matrix J has been destroyed.  
 21:    The operation MatMult(A,D2,D1) must be well defined.
 22:    Before calling the operation MatGetDiagonal(), the function 
 23:    MatADAComputeDiagonal() must be called.  The matrices A and D1 must
 24:    be the same during calls to MatADAComputeDiagonal() and
 25:    MatGetDiagonal().

 27:    Level: developer

 29: .seealso: MatCreate()
 30: @*/
 31: int MatCreateADA(Mat mat,Vec D1, Vec D2, Mat *J)
 32: {
 33:   MPI_Comm     comm=mat->comm;
 34:   TaoMatADACtx ctx;
 35:   int          info,nloc,n;

 38:   /*
 39:   info=MatCheckVecs(mat,D1,D2,&flg);CHKERRQ(info);
 40:   if (flg==PETSC_FALSE){
 41:     SETERRQ(PETSC_ERR_SUP,"InCompatible matrix and vector for ADA^T matrix");
 42:   }
 43:   */
 44:   info = PetscNew(_p_TaoMatADACtx,&ctx);CHKERRQ(info);

 46:   ctx->A=mat;
 47:   ctx->D1=D1;
 48:   ctx->D2=D2;
 49:   if (D1){
 50:     info = VecDuplicate(D1,&ctx->W);CHKERRQ(info);
 51:     info =  PetscObjectReference((PetscObject)D1);CHKERRQ(info);
 52:   } else {
 53:     ctx->W=0;
 54:   }
 55:   if (D2){
 56:     info = VecDuplicate(D2,&ctx->W2);CHKERRQ(info);
 57:     info = VecDuplicate(D2,&ctx->ADADiag);CHKERRQ(info);
 58:     info =  PetscObjectReference((PetscObject)D2);CHKERRQ(info);
 59:   } else {
 60:     ctx->W2=0;
 61:     ctx->ADADiag=0;
 62:   }

 64:   ctx->GotDiag=0;
 65:   info =  PetscObjectReference((PetscObject)mat);CHKERRQ(info);

 67:   info=VecGetLocalSize(D2,&nloc);CHKERRQ(info);
 68:   info=VecGetSize(D2,&n);CHKERRQ(info);

 70:   info = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(info);

 72:   info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_ADA);CHKERRQ(info);
 73:   info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_ADA);CHKERRQ(info);
 74:   info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_ADA);CHKERRQ(info);
 75:   info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_ADA);CHKERRQ(info);
 76:   info = MatShellSetOperation(*J,MATOP_DIAGONAL_SHIFT,(void(*)())MatDiagonalShift_ADA);CHKERRQ(info);
 77:   info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_ADA);CHKERRQ(info);
 78:   info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_ADA);CHKERRQ(info);
 79:   info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_ADA);CHKERRQ(info);
 80:   info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_ADA);CHKERRQ(info);
 81:   info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_ADA);CHKERRQ(info);
 82:   info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_ADA);CHKERRQ(info);
 83:   info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_ADA);CHKERRQ(info);
 84:   info = MatShellSetOperation(*J,MATOP_DUPLCIATE,(void(*)())MatDuplicate_ADA);CHKERRQ(info);
 85:   info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_ADA);CHKERRQ(info);

 87:   info = PetscLogObjectParent(*J,ctx->W); CHKERRQ(info);
 88:   info = PetscLogObjectParent(mat,*J); CHKERRQ(info);

 90:   info = MatSetOption(*J,MAT_SYMMETRIC);CHKERRQ(info);
 91:   return(0);
 92: }

 96: int MatMult_ADA(Mat mat,Vec a,Vec y)
 97: {
 98:   TaoMatADACtx ctx;
 99:   PetscScalar        one = 1.0;
100:   int           info;

103:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

105:   info = MatMult(ctx->A,a,ctx->W);CHKERRQ(info);
106:   if (ctx->D1){
107:     info = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(info);
108:   }
109:   info = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(info);
110:   if (ctx->D2){
111:     info = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(info);
112:     info = VecAXPY(y, one, ctx->W2);CHKERRQ(info);
113:   }
114:   return(0);
115: }

119: int MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
120: {
121:   int info;

124:   info = MatMult_ADA(mat,a,y);CHKERRQ(info);
125:   return(0);
126: }

130: int MatDiagonalShift_ADA(Vec D, Mat M)
131: {
132:   TaoMatADACtx ctx;
133:   PetscScalar        zero=0.0,one = 1.0;
134:   int           info;

137:   info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);

139:   if (ctx->D2==PETSC_NULL){
140:     info = VecDuplicate(D,&ctx->D2);CHKERRQ(info);
141:     info = VecSet(ctx->D2, zero);CHKERRQ(info);
142:   }
143:   info = VecAXPY(ctx->D2, one, D);CHKERRQ(info);

145:   return(0);
146: }

150: int MatDestroy_ADA(Mat mat)
151: {
152:   int          info;
153:   TaoMatADACtx ctx;

156:   info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
157:   info=VecDestroy(ctx->W);CHKERRQ(info);
158:   info=VecDestroy(ctx->W2);CHKERRQ(info);
159:   info=VecDestroy(ctx->ADADiag);CHKERRQ(info);
160:   info=MatDestroy(ctx->A);CHKERRQ(info);
161:   info=VecDestroy(ctx->D1);CHKERRQ(info);
162:   info=VecDestroy(ctx->D2);CHKERRQ(info);
163:   PetscFree(ctx);
164:   return(0);
165: }

169: int MatView_ADA(Mat mat,PetscViewer viewer)
170: {

173:   /*
174:   info = ViewerGetFormat(viewer,&format);CHKERRQ(info);
175:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) {
176:     return(0);  / * do nothing for now * /
177:   }
178:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
179:   info = MatView(ctx->A,viewer);CHKERRQ(info);
180:   if (ctx->D1){
181:     info = VecView(ctx->D1,viewer);CHKERRQ(info);
182:   }
183:   if (ctx->D2){
184:     info = VecView(ctx->D2,viewer);CHKERRQ(info);
185:   }
186:   */
187:   return(0);
188: }

192: int MatShift_ADA(Mat Y, PetscScalar a)
193: {
194:   int          info;
195:   TaoMatADACtx ctx;

198:   info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
199:   info = VecShift(ctx->D2,a);CHKERRQ(info);
200:   return(0);
201: }

205: int MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
206: {
207:   int          info;
208:   TaoMatADACtx ctx;
209:   Mat          A2;
210:   Vec          D1b=NULL,D2b;

213:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
214:   info = MatDuplicate(ctx->A,op,&A2);CHKERRQ(info);
215:   if (ctx->D1){
216:     info = VecDuplicate(ctx->D1,&D1b);CHKERRQ(info);
217:     info = VecCopy(ctx->D1,D1b);CHKERRQ(info);
218:   }
219:   info = VecDuplicate(ctx->D2,&D2b);CHKERRQ(info);
220:   info = VecCopy(ctx->D2,D2b);CHKERRQ(info);
221:   info = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(info);
222:   if (ctx->D1){
223:   info = PetscObjectDereference((PetscObject)D1b);CHKERRQ(info);
224:   }
225:   info = PetscObjectDereference((PetscObject)D2b);CHKERRQ(info);
226:   info = PetscObjectDereference((PetscObject)A2);CHKERRQ(info);

228:   return(0);
229: }

233: int MatEqual_ADA(Mat A,Mat B,PetscTruth *flg)
234: {
235:   int          info;
236:   TaoMatADACtx  ctx1,ctx2;

239:   info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
240:   info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
241:   info = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(info);
242:   if (*flg==PETSC_TRUE){
243:     info = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(info);
244:   }
245:   if (*flg==PETSC_TRUE){
246:     info = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(info);
247:   }
248:   return(0);
249: }

253: int MatScale_ADA(Mat mat, PetscScalar a)
254: {
255:   int          info;
256:   TaoMatADACtx ctx;

259:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
260:   info = VecScale(ctx->D1,a);CHKERRQ(info);
261:   if (ctx->D2){
262:     info = VecScale(ctx->D2,a);CHKERRQ(info);
263:   }
264:   return(0);
265: }

269: int MatTranspose_ADA(Mat mat,Mat *B)
270: {
271:   int          info;
272:   TaoMatADACtx ctx;

275:   if (*B){
276:     info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
277:     info = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(info);
278:   }
279:   return(0);
280: }

284: int MatADAComputeDiagonal(Mat mat)
285: {
286:   int          i,m,n,low,high,info;
287:   PetscScalar       *dtemp,*dptr;
288:   TaoMatADACtx ctx;

291:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

293:   info = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(info);
294:   info = MatGetSize(mat,&m,&n);CHKERRQ(info);
295: 
296:   info = PetscMalloc( n*sizeof(PetscScalar),&dtemp ); CHKERRQ(info);

298:   for (i=0; i<n; i++){
299:     info = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(info);
300:     info = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(info);
301:     info = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
302:   }
303:   for (i=0; i<n; i++){
304:     info = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
305:   }

307:   info = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(info);
308:   for (i=low; i<high; i++){
309:     dptr[i-low]= dtemp[i];
310:   }
311:   info = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(info);
312:   if (dtemp) PetscFree(dtemp);
313:   return(0);
314: }

318: int MatGetDiagonal_ADA(Mat mat,Vec v)
319: {
320:   int          info;
321:   PetscScalar       one=1.0;
322:   TaoMatADACtx ctx;

325:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
326:   info = MatADAComputeDiagonal(mat);
327:   info=VecCopy(ctx->ADADiag,v);CHKERRQ(info);
328:   if (ctx->D2){
329:     info=VecAXPY(v, one, ctx->D2);CHKERRQ(info);
330:   }

332:   return(0);
333: }

337: int MatGetSubMatrices_ADA(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
338: {
339:   int info,i;

342:   if (scall == MAT_INITIAL_MATRIX) {
343:     info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
344:   }

346:   for ( i=0; i<n; i++ ) {
347:     info = MatGetSubMatrix_ADA(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(info);
348:   }
349:   return(0);
350: }

354: int MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,
355:                         Mat *newmat)
356: {
357:   int          info,low,high;
358:   int          n,nlocal,i;
359:   int          *iptr;
360:   PetscScalar       *dptr,*ddptr,zero=0.0;
361:   VecType type_name;
362:   IS           ISrow;
363:   Vec          D1,D2;
364:   Mat          Atemp;
365:   TaoMatADACtx ctx;


369:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

371:   info = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(info);
372:   info = ISCreateStride(mat->comm,high-low,low,1,&ISrow);CHKERRQ(info);
373:   info = MatGetSubMatrix(ctx->A,ISrow,iscol,csize,cll,&Atemp);CHKERRQ(info);
374:   info = ISDestroy(ISrow);CHKERRQ(info);

376:   if (ctx->D1){
377:     info=VecDuplicate(ctx->D1,&D1);CHKERRQ(info);
378:     info=VecCopy(ctx->D1,D1);CHKERRQ(info);
379:   } else {
380:     D1=PETSC_NULL;
381:   }

383:   if (ctx->D2){
384:     info=ISGetSize(isrow,&n);CHKERRQ(info);
385:     info=ISGetLocalSize(isrow,&nlocal);CHKERRQ(info);
386:     info=VecCreate(ctx->D2->comm,&D2);CHKERRQ(info);
387:     info=VecGetType(ctx->D2,&type_name);CHKERRQ(info);
388:     info=VecSetSizes(D2,nlocal,n);CHKERRQ(info);
389:     info=VecSetType(D2,type_name);CHKERRQ(info);
390:     info=VecSet(D2, zero);CHKERRQ(info);
391:     info=VecGetArray(ctx->D2, &dptr); CHKERRQ(info);
392:     info=VecGetArray(D2, &ddptr); CHKERRQ(info);
393:     info=ISGetIndices(isrow,&iptr); CHKERRQ(info);
394:     for (i=0;i<nlocal;i++){
395:       ddptr[i] = dptr[iptr[i]-low];
396:     }
397:     info=ISRestoreIndices(isrow,&iptr); CHKERRQ(info);
398:     info=VecRestoreArray(D2, &ddptr); CHKERRQ(info);
399:     info=VecRestoreArray(ctx->D2, &dptr); CHKERRQ(info);
400: 
401:   } else {
402:     D2=PETSC_NULL;
403:   }

405:   info = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(info);
406:   info = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(info);
407:   info = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(info);
408:   if (ctx->D1){
409:     info = PetscObjectDereference((PetscObject)D1);CHKERRQ(info);
410:   }
411:   if (ctx->D2){
412:     info = PetscObjectDereference((PetscObject)D2);CHKERRQ(info);
413:   }
414:   return(0);
415: }

419: int MatGetRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
420: {
421:   int info,m,n;

424:   info = MatGetSize(mat,&m,&n);CHKERRQ(info);

426:   if (*ncols>0){
427:     info = PetscMalloc( (*ncols)*sizeof(int),cols );CHKERRQ(info);
428:     info = PetscMalloc( (*ncols)*sizeof(PetscScalar),vals );CHKERRQ(info);
429:   } else {
430:     *cols=PETSC_NULL;
431:     *vals=PETSC_NULL;
432:   }
433: 
434:   return(0);
435: }

439: int MatRestoreRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
440: {
442:   if (*ncols>0){
443:     PetscFree(*cols);
444:     PetscFree(*vals);
445:   }
446:   *cols=PETSC_NULL;
447:   *vals=PETSC_NULL;
448:   return(0);
449: }

453: int MatGetColumnVector_ADA(Mat mat,Vec Y, int col)
454: {
455:   int    info,low,high;
456:   PetscScalar zero=0.0,one=1.0;

459:   info=VecSet(Y, zero);CHKERRQ(info);
460:   info=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(info);
461:   if (col>=low && col<high){
462:     info=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(info);
463:   }
464:   info=VecAssemblyBegin(Y);CHKERRQ(info);
465:   info=VecAssemblyEnd(Y);CHKERRQ(info);
466:   info=MatMult_ADA(mat,Y,Y);CHKERRQ(info);

468:   return(0);
469: }

471: int MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
472: {
473:   int info,size;
474:   TaoMatADACtx  ctx;

477:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
478:   MPI_Comm_size(mat->comm,&size);

480:   if (newtype==MATSAME){

482:     info=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(info);

484:   } else if (newtype==MATMPIDENSE){

486:     int i,j,low,high,m,n,M,N;
487:     PetscScalar *dptr;
488:     Vec X;

490:     info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
491:     info=MatGetSize(mat,&M,&N);CHKERRQ(info);
492:     info=MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
493:     info = MatCreateMPIDense(mat->comm,m,m,N,N,PETSC_NULL,NewMat);
494:     CHKERRQ(info);
495:     info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
496:     for (i=0;i<M;i++){
497:       info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
498:       info = VecGetArray(X,&dptr);CHKERRQ(info);
499:       for (j=0; j<high-low; j++){
500:         info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
501:       }
502:       info=VecRestoreArray(X,&dptr);CHKERRQ(info);
503:     }
504:     info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
505:     info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
506:     info = VecDestroy(X);CHKERRQ(info);

508:   } else if (newtype==MATSEQDENSE && size==1){

510:     int i,j,low,high,m,n,M,N;
511:     PetscScalar *dptr;
512:     Vec X;

514:     info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
515:     info = MatGetSize(mat,&M,&N);CHKERRQ(info);
516:     info = MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
517:     info = MatCreateSeqDense(mat->comm,N,N,PETSC_NULL,NewMat);
518:     CHKERRQ(info);
519:     info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
520:     for (i=0;i<M;i++){
521:       info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
522:       info = VecGetArray(X,&dptr);CHKERRQ(info);
523:       for (j=0; j<high-low; j++){
524:         info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
525:       }
526:       info=VecRestoreArray(X,&dptr);CHKERRQ(info);
527:     }
528:     info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
529:     info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
530:     info=VecDestroy(X);CHKERRQ(info);

532:   } else {
533:     SETERRQ(1,"No support to convert objects to that type");
534:   }
535:   return(0);
536: }

540: int MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
541: {
542:   int info;
543:   TaoMatADACtx  ctx;

546:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

548:   if (type == NORM_FROBENIUS) {
549:     *norm = 1.0;
550:   } else if (type == NORM_1 || type == NORM_INFINITY) {
551:     *norm = 1.0;
552:   } else {
553:     SETERRQ(PETSC_ERR_SUP,"No two norm");
554:   }
555:   return(0);
556: }