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: }