Actual source code: axpy.c

  1: #define PETSCMAT_DLL

 3:  #include include/private/matimpl.h

  7: /*@
  8:    MatAXPY - Computes Y = a*X + Y.

 10:    Collective on Mat

 12:    Input Parameters:
 13: +  a - the scalar multiplier
 14: .  X - the first matrix
 15: .  Y - the second matrix
 16: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN

 18:    Notes:
 19:      Will only be efficient if one has the SAME_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN

 21:    Level: intermediate

 23: .keywords: matrix, add

 25: .seealso: MatAYPX()
 26:  @*/
 27: PetscErrorCode  MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)
 28: {
 30:   PetscInt       m1,m2,n1,n2;


 36:   MatGetSize(X,&m1,&n1);
 37:   MatGetSize(Y,&m2,&n2);
 38:   if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %D %D %D %D",m1,m2,n1,n2);

 40:   if (Y->ops->axpy) {
 41:     (*Y->ops->axpy)(Y,a,X,str);
 42:   } else {
 43:     MatAXPY_Basic(Y,a,X,str);
 44:   }
 45:   return(0);
 46: }


 51: PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
 52: {
 53:   PetscInt          i,start,end,j,ncols,m,n;
 54:   PetscErrorCode    ierr;
 55:   const PetscInt    *row;
 56:   PetscScalar       *val;
 57:   const PetscScalar *vals;

 60:   MatGetSize(X,&m,&n);
 61:   MatGetOwnershipRange(X,&start,&end);
 62:   if (a == 1.0) {
 63:     for (i = start; i < end; i++) {
 64:       MatGetRow(X,i,&ncols,&row,&vals);
 65:       MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);
 66:       MatRestoreRow(X,i,&ncols,&row,&vals);
 67:     }
 68:   } else {
 69:     PetscMalloc((n+1)*sizeof(PetscScalar),&val);
 70:     for (i=start; i<end; i++) {
 71:       MatGetRow(X,i,&ncols,&row,&vals);
 72:       for (j=0; j<ncols; j++) {
 73:         val[j] = a*vals[j];
 74:       }
 75:       MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);
 76:       MatRestoreRow(X,i,&ncols,&row,&vals);
 77:     }
 78:     PetscFree(val);
 79:   }
 80:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
 81:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
 82:   return(0);
 83: }

 87: /*@
 88:    MatShift - Computes Y =  Y + a I, where a is a PetscScalar and I is the identity matrix.

 90:    Collective on Mat

 92:    Input Parameters:
 93: +  Y - the matrices
 94: -  a - the PetscScalar 

 96:    Level: intermediate

 98: .keywords: matrix, add, shift

100: .seealso: MatDiagonalSet()
101:  @*/
102: PetscErrorCode  MatShift(Mat Y,PetscScalar a)
103: {
105:   PetscInt       i,start,end;

109:   if (!Y->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
110:   if (Y->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
111:   MatPreallocated(Y);

113:   if (Y->ops->shift) {
114:     (*Y->ops->shift)(Y,a);
115:   } else {
116:     PetscScalar alpha = a;
117:     MatGetOwnershipRange(Y,&start,&end);
118:     for (i=start; i<end; i++) {
119:       MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
120:     }
121:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
122:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
123:   }
124:   return(0);
125: }

129: PetscErrorCode  MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is)
130: {
132:   PetscInt       i,start,end,vstart,vend;
133:   PetscScalar    *v;

136:   VecGetOwnershipRange(D,&vstart,&vend);
137:   MatGetOwnershipRange(Y,&start,&end);
138:   if (vstart != start || vend != end) {
139:     SETERRQ4(PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %D %D vec %D %D mat",vstart,vend,start,end);
140:   }
141:   VecGetArray(D,&v);
142:   for (i=start; i<end; i++) {
143:     MatSetValues(Y,1,&i,1,&i,v+i-start,is);
144:   }
145:   VecRestoreArray(D,&v);
146:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
147:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
148:   return(0);
149: }

153: /*@
154:    MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
155:    that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
156:    INSERT_VALUES.

158:    Input Parameters:
159: +  Y - the input matrix
160: .  D - the diagonal matrix, represented as a vector
161: -  i - INSERT_VALUES or ADD_VALUES

163:    Collective on Mat and Vec

165:    Level: intermediate

167: .keywords: matrix, add, shift, diagonal

169: .seealso: MatShift()
170: @*/
171: PetscErrorCode  MatDiagonalSet(Mat Y,Vec D,InsertMode is)
172: {

178:   if (Y->ops->diagonalset) {
179:     (*Y->ops->diagonalset)(Y,D,is);
180:   } else {
181:     MatDiagonalSet_Default(Y,D,is);
182:   }
183:   return(0);
184: }

188: /*@
189:    MatAYPX - Computes Y = a*Y + X.

191:    Collective on Mat

193:    Input Parameters:
194: +  a - the PetscScalar multiplier
195: .  Y - the first matrix
196: .  X - the second matrix
197: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN 

199:    Notes:
200:      Will only be efficient if one has the SAME_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN

202:    Level: intermediate

204: .keywords: matrix, add

206: .seealso: MatAXPY()
207:  @*/
208: PetscErrorCode  MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str)
209: {
210:   PetscScalar    one = 1.0;
212:   PetscInt       mX,mY,nX,nY;


218:   MatGetSize(X,&mX,&nX);
219:   MatGetSize(X,&mY,&nY);
220:   if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrices: %D %D first %D %D second",mX,mY,nX,nY);

222:   MatScale(Y,a);
223:   MatAXPY(Y,one,X,str);CHKERRQ(ierr)
224:   return(0);
225: }

229: /*@
230:     MatComputeExplicitOperator - Computes the explicit matrix

232:     Collective on Mat

234:     Input Parameter:
235: .   inmat - the matrix

237:     Output Parameter:
238: .   mat - the explict preconditioned operator

240:     Notes:
241:     This computation is done by applying the operators to columns of the 
242:     identity matrix.

244:     Currently, this routine uses a dense matrix format when 1 processor
245:     is used and a sparse format otherwise.  This routine is costly in general,
246:     and is recommended for use only with relatively small systems.

248:     Level: advanced
249:    
250: .keywords: Mat, compute, explicit, operator

252: @*/
253: PetscErrorCode  MatComputeExplicitOperator(Mat inmat,Mat *mat)
254: {
255:   Vec            in,out;
257:   PetscInt       i,M,m,*rows,start,end;
258:   MPI_Comm       comm;
259:   PetscScalar    *array,zero = 0.0,one = 1.0;
260:   PetscMPIInt    size;


266:   comm = inmat->comm;
267:   MPI_Comm_size(comm,&size);

269:   MatGetLocalSize(inmat,&m,0);
270:   MatGetSize(inmat,&M,0);
271:   VecCreateMPI(comm,m,M,&in);
272:   VecDuplicate(in,&out);
273:   VecGetOwnershipRange(in,&start,&end);
274:   PetscMalloc((m+1)*sizeof(PetscInt),&rows);
275:   for (i=0; i<m; i++) {rows[i] = start + i;}

277:   MatCreate(comm,mat);
278:   MatSetSizes(*mat,m,m,M,M);
279:   if (size == 1) {
280:     MatSetType(*mat,MATSEQDENSE);
281:     MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
282:   } else {
283:     MatSetType(*mat,MATMPIAIJ);
284:     MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
285:   }

287:   for (i=0; i<M; i++) {

289:     VecSet(in,zero);
290:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
291:     VecAssemblyBegin(in);
292:     VecAssemblyEnd(in);

294:     MatMult(inmat,in,out);
295: 
296:     VecGetArray(out,&array);
297:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
298:     VecRestoreArray(out,&array);

300:   }
301:   PetscFree(rows);
302:   VecDestroy(out);
303:   VecDestroy(in);
304:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
305:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
306:   return(0);
307: }

309: /* Get the map xtoy which is used by MatAXPY() in the case of SUBSET_NONZERO_PATTERN */
312: PetscErrorCode MatAXPYGetxtoy_Private(PetscInt m,PetscInt *xi,PetscInt *xj,PetscInt *xgarray, PetscInt *yi,PetscInt *yj,PetscInt *ygarray, PetscInt **xtoy)
313: {
315:   PetscInt       row,i,nz,xcol,ycol,jx,jy,*x2y;

318:   PetscMalloc(xi[m]*sizeof(PetscInt),&x2y);
319:   i = 0;
320:   for (row=0; row<m; row++){
321:     nz = xi[1] - xi[0];
322:     jy = 0;
323:     for (jx=0; jx<nz; jx++,jy++){
324:       if (xgarray && ygarray){
325:         xcol = xgarray[xj[*xi + jx]];
326:         ycol = ygarray[yj[*yi + jy]];
327:       } else {
328:         xcol = xj[*xi + jx];
329:         ycol = yj[*yi + jy];  /* col index for y */
330:       }
331:       while ( ycol < xcol ) {
332:         jy++;
333:         if (ygarray){
334:           ycol = ygarray[yj[*yi + jy]];
335:         } else {
336:           ycol = yj[*yi + jy];
337:         }
338:       }
339:       if (xcol != ycol) SETERRQ2(PETSC_ERR_ARG_WRONG,"X matrix entry (%D,%D) is not in Y matrix",row,ycol);
340:       x2y[i++] = *yi + jy;
341:     }
342:     xi++; yi++;
343:   }
344:   *xtoy = x2y;
345:   return(0);
346: }