Actual source code: normm.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  2: #include <petsc/private/matimpl.h>          /*I "petscmat.h" I*/

  4: typedef struct {
  5:   Mat         A;
  6:   Vec         w,left,right,leftwork,rightwork;
  7:   PetscScalar scale;
  8: } Mat_Normal;

 12: PetscErrorCode MatScale_Normal(Mat inA,PetscScalar scale)
 13: {
 14:   Mat_Normal *a = (Mat_Normal*)inA->data;

 17:   a->scale *= scale;
 18:   return(0);
 19: }

 23: PetscErrorCode MatDiagonalScale_Normal(Mat inA,Vec left,Vec right)
 24: {
 25:   Mat_Normal     *a = (Mat_Normal*)inA->data;

 29:   if (left) {
 30:     if (!a->left) {
 31:       VecDuplicate(left,&a->left);
 32:       VecCopy(left,a->left);
 33:     } else {
 34:       VecPointwiseMult(a->left,left,a->left);
 35:     }
 36:   }
 37:   if (right) {
 38:     if (!a->right) {
 39:       VecDuplicate(right,&a->right);
 40:       VecCopy(right,a->right);
 41:     } else {
 42:       VecPointwiseMult(a->right,right,a->right);
 43:     }
 44:   }
 45:   return(0);
 46: }

 50: PetscErrorCode MatMult_Normal(Mat N,Vec x,Vec y)
 51: {
 52:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 54:   Vec            in;

 57:   in = x;
 58:   if (Na->right) {
 59:     if (!Na->rightwork) {
 60:       VecDuplicate(Na->right,&Na->rightwork);
 61:     }
 62:     VecPointwiseMult(Na->rightwork,Na->right,in);
 63:     in   = Na->rightwork;
 64:   }
 65:   MatMult(Na->A,in,Na->w);
 66:   MatMultTranspose(Na->A,Na->w,y);
 67:   if (Na->left) {
 68:     VecPointwiseMult(y,Na->left,y);
 69:   }
 70:   VecScale(y,Na->scale);
 71:   return(0);
 72: }

 76: PetscErrorCode MatMultAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
 77: {
 78:   Mat_Normal     *Na = (Mat_Normal*)N->data;
 80:   Vec            in;

 83:   in = v1;
 84:   if (Na->right) {
 85:     if (!Na->rightwork) {
 86:       VecDuplicate(Na->right,&Na->rightwork);
 87:     }
 88:     VecPointwiseMult(Na->rightwork,Na->right,in);
 89:     in   = Na->rightwork;
 90:   }
 91:   MatMult(Na->A,in,Na->w);
 92:   VecScale(Na->w,Na->scale);
 93:   if (Na->left) {
 94:     MatMultTranspose(Na->A,Na->w,v3);
 95:     VecPointwiseMult(v3,Na->left,v3);
 96:     VecAXPY(v3,1.0,v2);
 97:   } else {
 98:     MatMultTransposeAdd(Na->A,Na->w,v2,v3);
 99:   }
100:   return(0);
101: }

105: PetscErrorCode MatMultTranspose_Normal(Mat N,Vec x,Vec y)
106: {
107:   Mat_Normal     *Na = (Mat_Normal*)N->data;
109:   Vec            in;

112:   in = x;
113:   if (Na->left) {
114:     if (!Na->leftwork) {
115:       VecDuplicate(Na->left,&Na->leftwork);
116:     }
117:     VecPointwiseMult(Na->leftwork,Na->left,in);
118:     in   = Na->leftwork;
119:   }
120:   MatMult(Na->A,in,Na->w);
121:   MatMultTranspose(Na->A,Na->w,y);
122:   if (Na->right) {
123:     VecPointwiseMult(y,Na->right,y);
124:   }
125:   VecScale(y,Na->scale);
126:   return(0);
127: }

131: PetscErrorCode MatMultTransposeAdd_Normal(Mat N,Vec v1,Vec v2,Vec v3)
132: {
133:   Mat_Normal     *Na = (Mat_Normal*)N->data;
135:   Vec            in;

138:   in = v1;
139:   if (Na->left) {
140:     if (!Na->leftwork) {
141:       VecDuplicate(Na->left,&Na->leftwork);
142:     }
143:     VecPointwiseMult(Na->leftwork,Na->left,in);
144:     in   = Na->leftwork;
145:   }
146:   MatMult(Na->A,in,Na->w);
147:   VecScale(Na->w,Na->scale);
148:   if (Na->right) {
149:     MatMultTranspose(Na->A,Na->w,v3);
150:     VecPointwiseMult(v3,Na->right,v3);
151:     VecAXPY(v3,1.0,v2);
152:   } else {
153:     MatMultTransposeAdd(Na->A,Na->w,v2,v3);
154:   }
155:   return(0);
156: }

160: PetscErrorCode MatDestroy_Normal(Mat N)
161: {
162:   Mat_Normal     *Na = (Mat_Normal*)N->data;

166:   MatDestroy(&Na->A);
167:   VecDestroy(&Na->w);
168:   VecDestroy(&Na->left);
169:   VecDestroy(&Na->right);
170:   VecDestroy(&Na->leftwork);
171:   VecDestroy(&Na->rightwork);
172:   PetscFree(N->data);
173:   return(0);
174: }

176: /*
177:       Slow, nonscalable version
178: */
181: PetscErrorCode MatGetDiagonal_Normal(Mat N,Vec v)
182: {
183:   Mat_Normal        *Na = (Mat_Normal*)N->data;
184:   Mat               A   = Na->A;
185:   PetscErrorCode    ierr;
186:   PetscInt          i,j,rstart,rend,nnz;
187:   const PetscInt    *cols;
188:   PetscScalar       *diag,*work,*values;
189:   const PetscScalar *mvalues;

192:   PetscMalloc2(A->cmap->N,&diag,A->cmap->N,&work);
193:   PetscMemzero(work,A->cmap->N*sizeof(PetscScalar));
194:   MatGetOwnershipRange(A,&rstart,&rend);
195:   for (i=rstart; i<rend; i++) {
196:     MatGetRow(A,i,&nnz,&cols,&mvalues);
197:     for (j=0; j<nnz; j++) {
198:       work[cols[j]] += mvalues[j]*mvalues[j];
199:     }
200:     MatRestoreRow(A,i,&nnz,&cols,&mvalues);
201:   }
202:   MPI_Allreduce(work,diag,A->cmap->N,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
203:   rstart = N->cmap->rstart;
204:   rend   = N->cmap->rend;
205:   VecGetArray(v,&values);
206:   PetscMemcpy(values,diag+rstart,(rend-rstart)*sizeof(PetscScalar));
207:   VecRestoreArray(v,&values);
208:   PetscFree2(diag,work);
209:   VecScale(v,Na->scale);
210:   return(0);
211: }

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

218:    Collective on Mat

220:    Input Parameter:
221: .   A  - the (possibly rectangular) matrix

223:    Output Parameter:
224: .   N - the matrix that represents A'*A

226:    Level: intermediate

228:    Notes: The product A'*A is NOT actually formed! Rather the new matrix
229:           object performs the matrix-vector product by first multiplying by
230:           A and then A'
231: @*/
232: PetscErrorCode  MatCreateNormal(Mat A,Mat *N)
233: {
235:   PetscInt       m,n;
236:   Mat_Normal     *Na;

239:   MatGetLocalSize(A,&m,&n);
240:   MatCreate(PetscObjectComm((PetscObject)A),N);
241:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
242:   PetscObjectChangeTypeName((PetscObject)*N,MATNORMAL);

244:   PetscNewLog(*N,&Na);
245:   (*N)->data = (void*) Na;
246:   PetscObjectReference((PetscObject)A);
247:   Na->A      = A;
248:   Na->scale  = 1.0;

250:   VecCreateMPI(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,&Na->w);

252:   (*N)->ops->destroy          = MatDestroy_Normal;
253:   (*N)->ops->mult             = MatMult_Normal;
254:   (*N)->ops->multtranspose    = MatMultTranspose_Normal;
255:   (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal;
256:   (*N)->ops->multadd          = MatMultAdd_Normal;
257:   (*N)->ops->getdiagonal      = MatGetDiagonal_Normal;
258:   (*N)->ops->scale            = MatScale_Normal;
259:   (*N)->ops->diagonalscale    = MatDiagonalScale_Normal;
260:   (*N)->assembled             = PETSC_TRUE;
261:   (*N)->cmap->N               = A->cmap->N;
262:   (*N)->rmap->N               = A->cmap->N;
263:   (*N)->cmap->n               = A->cmap->n;
264:   (*N)->rmap->n               = A->cmap->n;
265:   return(0);
266: }