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;
14: a->scale *= scale;
15: return 0;
16: }
18: PetscErrorCode MatDiagonalScale_Normal(Mat inA,Vec left,Vec right)
19: {
20: Mat_Normal *a = (Mat_Normal*)inA->data;
22: if (left) {
23: if (!a->left) {
24: VecDuplicate(left,&a->left);
25: VecCopy(left,a->left);
26: } else {
27: VecPointwiseMult(a->left,left,a->left);
28: }
29: }
30: if (right) {
31: if (!a->right) {
32: VecDuplicate(right,&a->right);
33: VecCopy(right,a->right);
34: } else {
35: VecPointwiseMult(a->right,right,a->right);
36: }
37: }
38: return 0;
39: }
41: PetscErrorCode MatIncreaseOverlap_Normal(Mat A,PetscInt is_max,IS is[],PetscInt ov)
42: {
43: Mat_Normal *a = (Mat_Normal*)A->data;
44: Mat pattern;
47: MatProductCreate(a->A,a->A,NULL,&pattern);
48: MatProductSetType(pattern,MATPRODUCT_AtB);
49: MatProductSetFromOptions(pattern);
50: MatProductSymbolic(pattern);
51: MatIncreaseOverlap(pattern,is_max,is,ov);
52: MatDestroy(&pattern);
53: return 0;
54: }
56: PetscErrorCode MatCreateSubMatrices_Normal(Mat mat,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[])
57: {
58: Mat_Normal *a = (Mat_Normal*)mat->data;
59: Mat B = a->A, *suba;
60: IS *row;
61: PetscInt M;
64: if (scall != MAT_REUSE_MATRIX) {
65: PetscCalloc1(n,submat);
66: }
67: MatGetSize(B,&M,NULL);
68: PetscMalloc1(n,&row);
69: ISCreateStride(PETSC_COMM_SELF,M,0,1,&row[0]);
70: ISSetIdentity(row[0]);
71: for (M = 1; M < n; ++M) row[M] = row[0];
72: MatCreateSubMatrices(B,n,row,icol,MAT_INITIAL_MATRIX,&suba);
73: for (M = 0; M < n; ++M) {
74: MatCreateNormal(suba[M],*submat+M);
75: ((Mat_Normal*)(*submat)[M]->data)->scale = a->scale;
76: }
77: ISDestroy(&row[0]);
78: PetscFree(row);
79: MatDestroySubMatrices(n,&suba);
80: return 0;
81: }
83: PetscErrorCode MatPermute_Normal(Mat A,IS rowp,IS colp,Mat *B)
84: {
85: Mat_Normal *a = (Mat_Normal*)A->data;
86: Mat C,Aa = a->A;
87: IS row;
90: ISCreateStride(PetscObjectComm((PetscObject)Aa),Aa->rmap->n,Aa->rmap->rstart,1,&row);
91: ISSetIdentity(row);
92: MatPermute(Aa,row,colp,&C);
93: ISDestroy(&row);
94: MatCreateNormal(C,B);
95: MatDestroy(&C);
96: return 0;
97: }
99: PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
100: {
101: Mat_Normal *a = (Mat_Normal*)A->data;
102: Mat C;
105: MatDuplicate(a->A,op,&C);
106: MatCreateNormal(C,B);
107: MatDestroy(&C);
108: if (op == MAT_COPY_VALUES) ((Mat_Normal*)(*B)->data)->scale = a->scale;
109: return 0;
110: }
112: PetscErrorCode MatCopy_Normal(Mat A,Mat B,MatStructure str)
113: {
114: Mat_Normal *a = (Mat_Normal*)A->data,*b = (Mat_Normal*)B->data;
117: MatCopy(a->A,b->A,str);
118: b->scale = a->scale;
119: VecDestroy(&b->left);
120: VecDestroy(&b->right);
121: VecDestroy(&b->leftwork);
122: VecDestroy(&b->rightwork);
123: return 0;
124: }
126: PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
127: {
128: Mat_Normal *Na = (Mat_Normal*)N->data;
129: Vec in;
131: in = x;
132: if (Na->right) {
133: if (!Na->rightwork) {
134: VecDuplicate(Na->right,&Na->rightwork);
135: }
136: VecPointwiseMult(Na->rightwork,Na->right,in);
137: in = Na->rightwork;
138: }
139: MatMult(Na->A,in,Na->w);
140: MatMultTranspose(Na->A,Na->w,y);
141: if (Na->left) {
142: VecPointwiseMult(y,Na->left,y);
143: }
144: VecScale(y,Na->scale);
145: return 0;
146: }
148: PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
149: {
150: Mat_Normal *Na = (Mat_Normal*)N->data;
151: Vec in;
153: in = v1;
154: if (Na->right) {
155: if (!Na->rightwork) {
156: VecDuplicate(Na->right,&Na->rightwork);
157: }
158: VecPointwiseMult(Na->rightwork,Na->right,in);
159: in = Na->rightwork;
160: }
161: MatMult(Na->A,in,Na->w);
162: VecScale(Na->w,Na->scale);
163: if (Na->left) {
164: MatMultTranspose(Na->A,Na->w,v3);
165: VecPointwiseMult(v3,Na->left,v3);
166: VecAXPY(v3,1.0,v2);
167: } else {
168: MatMultTransposeAdd(Na->A,Na->w,v2,v3);
169: }
170: return 0;
171: }
173: PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y)
174: {
175: Mat_Normal *Na = (Mat_Normal*)N->data;
176: Vec in;
178: in = x;
179: if (Na->left) {
180: if (!Na->leftwork) {
181: VecDuplicate(Na->left,&Na->leftwork);
182: }
183: VecPointwiseMult(Na->leftwork,Na->left,in);
184: in = Na->leftwork;
185: }
186: MatMult(Na->A,in,Na->w);
187: MatMultTranspose(Na->A,Na->w,y);
188: if (Na->right) {
189: VecPointwiseMult(y,Na->right,y);
190: }
191: VecScale(y,Na->scale);
192: return 0;
193: }
195: PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
196: {
197: Mat_Normal *Na = (Mat_Normal*)N->data;
198: Vec in;
200: in = v1;
201: if (Na->left) {
202: if (!Na->leftwork) {
203: VecDuplicate(Na->left,&Na->leftwork);
204: }
205: VecPointwiseMult(Na->leftwork,Na->left,in);
206: in = Na->leftwork;
207: }
208: MatMult(Na->A,in,Na->w);
209: VecScale(Na->w,Na->scale);
210: if (Na->right) {
211: MatMultTranspose(Na->A,Na->w,v3);
212: VecPointwiseMult(v3,Na->right,v3);
213: VecAXPY(v3,1.0,v2);
214: } else {
215: MatMultTransposeAdd(Na->A,Na->w,v2,v3);
216: }
217: return 0;
218: }
220: PetscErrorCode MatDestroy_Normal(Mat N)
221: {
222: Mat_Normal *Na = (Mat_Normal*)N->data;
224: MatDestroy(&Na->A);
225: VecDestroy(&Na->w);
226: VecDestroy(&Na->left);
227: VecDestroy(&Na->right);
228: VecDestroy(&Na->leftwork);
229: VecDestroy(&Na->rightwork);
230: PetscFree(N->data);
231: PetscObjectComposeFunction((PetscObject)N,"MatNormalGetMat_C",NULL);
232: PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_seqaij_C",NULL);
233: PetscObjectComposeFunction((PetscObject)N,"MatConvert_normal_mpiaij_C",NULL);
234: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_seqdense_C",NULL);
235: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_mpidense_C",NULL);
236: PetscObjectComposeFunction((PetscObject)N,"MatProductSetFromOptions_normal_dense_C",NULL);
237: return 0;
238: }
240: /*
241: Slow, nonscalable version
242: */
243: PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
244: {
245: Mat_Normal *Na = (Mat_Normal*)N->data;
246: Mat A = Na->A;
247: PetscInt i,j,rstart,rend,nnz;
248: const PetscInt *cols;
249: PetscScalar *diag,*work,*values;
250: const PetscScalar *mvalues;
252: PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
253: PetscArrayzero(work,A->cmap->N);
254: MatGetOwnershipRange(A,&rstart,&rend);
255: for (i=rstart; i<rend; i++) {
256: MatGetRow(A,i,&nnz,&cols,&mvalues);
257: for (j=0; j<nnz; j++) {
258: work[cols[j]] += mvalues[j]*mvalues[j];
259: }
260: MatRestoreRow(A,i,&nnz,&cols,&mvalues);
261: }
262: MPIU_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
263: rstart = N->cmap->rstart;
264: rend = N->cmap->rend;
265: VecGetArray(v,&values);
266: PetscArraycpy(values,diag+rstart,rend-rstart);
267: VecRestoreArray(v,&values);
268: PetscFree2(diag,work);
269: VecScale(v,Na->scale);
270: return 0;
271: }
273: PetscErrorCode MatNormalGetMat_Normal(Mat A,Mat *M)
274: {
275: Mat_Normal *Aa = (Mat_Normal*)A->data;
277: *M = Aa->A;
278: return 0;
279: }
281: /*@
282: MatNormalGetMat - Gets the Mat object stored inside a MATNORMAL
284: Logically collective on Mat
286: Input Parameter:
287: . A - the MATNORMAL matrix
289: Output Parameter:
290: . M - the matrix object stored inside A
292: Level: intermediate
294: .seealso: MatCreateNormal()
296: @*/
297: PetscErrorCode MatNormalGetMat(Mat A,Mat *M)
298: {
302: PetscUseMethod(A,"MatNormalGetMat_C",(Mat,Mat*),(A,M));
303: return 0;
304: }
306: PetscErrorCode MatConvert_Normal_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
307: {
308: Mat_Normal *Aa = (Mat_Normal*)A->data;
309: Mat B;
310: PetscInt m,n,M,N;
312: MatGetSize(A,&M,&N);
313: MatGetLocalSize(A,&m,&n);
314: if (reuse == MAT_REUSE_MATRIX) {
315: B = *newmat;
316: MatProductReplaceMats(Aa->A,Aa->A,NULL,B);
317: } else {
318: MatProductCreate(Aa->A,Aa->A,NULL,&B);
319: MatProductSetType(B,MATPRODUCT_AtB);
320: MatProductSetFromOptions(B);
321: MatProductSymbolic(B);
322: MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
323: }
324: MatProductNumeric(B);
325: if (reuse == MAT_INPLACE_MATRIX) {
326: MatHeaderReplace(A,&B);
327: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
328: MatConvert(*newmat,MATAIJ,MAT_INPLACE_MATRIX,newmat);
329: return 0;
330: }
332: typedef struct {
333: Mat work[2];
334: } Normal_Dense;
336: PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
337: {
338: Mat A,B;
339: Normal_Dense *contents;
340: Mat_Normal *a;
341: PetscScalar *array;
343: MatCheckProduct(C,3);
344: A = C->product->A;
345: a = (Mat_Normal*)A->data;
346: B = C->product->B;
347: contents = (Normal_Dense*)C->product->data;
349: if (a->right) {
350: MatCopy(B,C,SAME_NONZERO_PATTERN);
351: MatDiagonalScale(C,a->right,NULL);
352: }
353: MatProductNumeric(contents->work[0]);
354: MatDenseGetArrayWrite(C,&array);
355: MatDensePlaceArray(contents->work[1],array);
356: MatProductNumeric(contents->work[1]);
357: MatDenseRestoreArrayWrite(C,&array);
358: MatDenseResetArray(contents->work[1]);
359: MatSetOption(C,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
360: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
361: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
362: MatScale(C,a->scale);
363: return 0;
364: }
366: PetscErrorCode MatNormal_DenseDestroy(void *ctx)
367: {
368: Normal_Dense *contents = (Normal_Dense*)ctx;
370: MatDestroy(contents->work);
371: MatDestroy(contents->work+1);
372: PetscFree(contents);
373: return 0;
374: }
376: PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
377: {
378: Mat A,B;
379: Normal_Dense *contents = NULL;
380: Mat_Normal *a;
381: PetscScalar *array;
382: PetscInt n,N,m,M;
384: MatCheckProduct(C,4);
386: A = C->product->A;
387: a = (Mat_Normal*)A->data;
389: B = C->product->B;
390: MatGetLocalSize(C,&m,&n);
391: MatGetSize(C,&M,&N);
392: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
393: MatGetLocalSize(B,NULL,&n);
394: MatGetSize(B,NULL,&N);
395: MatGetLocalSize(A,&m,NULL);
396: MatGetSize(A,&M,NULL);
397: MatSetSizes(C,m,n,M,N);
398: }
399: MatSetType(C,((PetscObject)B)->type_name);
400: MatSetUp(C);
401: PetscNew(&contents);
402: C->product->data = contents;
403: C->product->destroy = MatNormal_DenseDestroy;
404: if (a->right) {
405: MatProductCreate(a->A,C,NULL,contents->work);
406: } else {
407: MatProductCreate(a->A,B,NULL,contents->work);
408: }
409: MatProductSetType(contents->work[0],MATPRODUCT_AB);
410: MatProductSetFromOptions(contents->work[0]);
411: MatProductSymbolic(contents->work[0]);
412: MatProductCreate(a->A,contents->work[0],NULL,contents->work+1);
413: MatProductSetType(contents->work[1],MATPRODUCT_AtB);
414: MatProductSetFromOptions(contents->work[1]);
415: MatProductSymbolic(contents->work[1]);
416: MatDenseGetArrayWrite(C,&array);
417: MatSeqDenseSetPreallocation(contents->work[1],array);
418: MatMPIDenseSetPreallocation(contents->work[1],array);
419: MatDenseRestoreArrayWrite(C,&array);
420: C->ops->productnumeric = MatProductNumeric_Normal_Dense;
421: return 0;
422: }
424: PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
425: {
426: C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
427: return 0;
428: }
430: PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
431: {
432: Mat_Product *product = C->product;
434: if (product->type == MATPRODUCT_AB) {
435: MatProductSetFromOptions_Normal_Dense_AB(C);
436: }
437: return 0;
438: }
440: /*@
441: MatCreateNormal - Creates a new matrix object that behaves like A'*A.
443: Collective on Mat
445: Input Parameter:
446: . A - the (possibly rectangular) matrix
448: Output Parameter:
449: . N - the matrix that represents A'*A
451: Level: intermediate
453: Notes:
454: The product A'*A is NOT actually formed! Rather the new matrix
455: object performs the matrix-vector product by first multiplying by
456: A and then A'
457: @*/
458: PetscErrorCode MatCreateNormal(Mat A,Mat *N)
459: {
460: PetscInt n,nn;
461: Mat_Normal *Na;
462: VecType vtype;
464: MatGetSize(A,NULL,&nn);
465: MatGetLocalSize(A,NULL,&n);
466: MatCreate(PetscObjectComm((PetscObject)A),N);
467: MatSetSizes(*N,n,n,nn,nn);
468: PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);
469: PetscLayoutReference(A->cmap,&(*N)->rmap);
470: PetscLayoutReference(A->cmap,&(*N)->cmap);
472: PetscNewLog(*N,&Na);
473: (*N)->data = (void*) Na;
474: PetscObjectReference((PetscObject)A);
475: Na->A = A;
476: Na->scale = 1.0;
478: MatCreateVecs(A,NULL,&Na->w);
480: (*N)->ops->destroy = MatDestroy_Normal;
481: (*N)->ops->mult = MatMult_Normal;
482: (*N)->ops->multtranspose = MatMultTranspose_Normal;
483: (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal;
484: (*N)->ops->multadd = MatMultAdd_Normal;
485: (*N)->ops->getdiagonal = MatGetDiagonal_Normal;
486: (*N)->ops->scale = MatScale_Normal;
487: (*N)->ops->diagonalscale = MatDiagonalScale_Normal;
488: (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal;
489: (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
490: (*N)->ops->permute = MatPermute_Normal;
491: (*N)->ops->duplicate = MatDuplicate_Normal;
492: (*N)->ops->copy = MatCopy_Normal;
493: (*N)->assembled = PETSC_TRUE;
494: (*N)->preallocated = PETSC_TRUE;
496: PetscObjectComposeFunction((PetscObject)(*N),"MatNormalGetMat_C",MatNormalGetMat_Normal);
497: PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_seqaij_C",MatConvert_Normal_AIJ);
498: PetscObjectComposeFunction((PetscObject)(*N),"MatConvert_normal_mpiaij_C",MatConvert_Normal_AIJ);
499: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_seqdense_C",MatProductSetFromOptions_Normal_Dense);
500: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_mpidense_C",MatProductSetFromOptions_Normal_Dense);
501: PetscObjectComposeFunction((PetscObject)(*N),"MatProductSetFromOptions_normal_dense_C",MatProductSetFromOptions_Normal_Dense);
502: MatSetOption(*N,MAT_SYMMETRIC,PETSC_TRUE);
503: MatGetVecType(A,&vtype);
504: MatSetVecType(*N,vtype);
505: #if defined(PETSC_HAVE_DEVICE)
506: MatBindToCPU(*N,A->boundtocpu);
507: #endif
508: return 0;
509: }