Actual source code: normm.c
petsc-3.8.4 2018-03-24
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: PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));
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: PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));
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: The product A'*A is NOT actually formed! Rather the new matrix
211: object performs the matrix-vector product by first multiplying by
212: A and then A'
213: @*/
214: PetscErrorCode MatCreateNormal(Mat A,Mat *N)
215: {
217: PetscInt m,n;
218: Mat_Normal *Na;
221: MatGetLocalSize(A,&m,&n);
222: MatCreate(PetscObjectComm((PetscObject)A),N);
223: MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
224: PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);
226: PetscNewLog(*N,&Na);
227: (*N)->data = (void*) Na;
228: PetscObjectReference((PetscObject)A);
229: Na->A = A;
230: Na->scale = 1.0;
232: VecCreateMPI(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,&Na->w);
234: (*N)->ops->destroy = MatDestroy_Normal;
235: (*N)->ops->mult = MatMult_Normal;
236: (*N)->ops->multtranspose = MatMultTranspose_Normal;
237: (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal;
238: (*N)->ops->multadd = MatMultAdd_Normal;
239: (*N)->ops->getdiagonal = MatGetDiagonal_Normal;
240: (*N)->ops->scale = MatScale_Normal;
241: (*N)->ops->diagonalscale = MatDiagonalScale_Normal;
242: (*N)->assembled = PETSC_TRUE;
243: (*N)->cmap->N = A->cmap->N;
244: (*N)->rmap->N = A->cmap->N;
245: (*N)->cmap->n = A->cmap->n;
246: (*N)->rmap->n = A->cmap->n;
248: (*N)->preallocated = PETSC_TRUE;
249: return(0);
250: }