1: /*$Id: axpy.c,v 1.54 2001/08/06 21:16:10 bsmith Exp bsmith $*/

 3:  #include src/mat/matimpl.h

  5: #undef __FUNCT__  
  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 or DIFFERENT_NONZERO_PATTERN

 18:    Contributed by: Matthew Knepley

 20:    Notes:
 21:      Will only be efficient if one has the SAME_NONZERO_PATTERN

 23:    Level: intermediate

 25: .keywords: matrix, add

 27: .seealso: MatAYPX()
 28:  @*/
 29: int MatAXPY(PetscScalar *a,Mat X,Mat Y,MatStructure str)
 30: {
 31:   int         m1,m2,n1,n2,ierr;


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

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


 51: #undef __FUNCT__  
 53: int MatAXPY_Basic(PetscScalar *a,Mat X,Mat Y,MatStructure str)
 54: {
 55:   int         i,*row,start,end,j,ncols,ierr,m,n;
 56:   PetscScalar *val,*vals;

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

 84: #undef __FUNCT__  
 86: /*@
 87:    MatShift - Computes Y =  Y + a I, where a is a PetscScalar and I is the identity matrix.

 89:    Collective on Mat

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

 95:    Level: intermediate

 97: .keywords: matrix, add, shift

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

108:   if (Y->ops->shift) {
109:     (*Y->ops->shift)(a,Y);
110:   } else {
111:     MatGetOwnershipRange(Y,&start,&end);
112:     for (i=start; i<end; i++) {
113:       MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES);
114:     }
115:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
116:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
117:   }
118:   return(0);
119: }

121: #undef __FUNCT__  
123: /*@
124:    MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
125:    that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
126:    INSERT_VALUES.

128:    Input Parameters:
129: +  Y - the input matrix
130: .  D - the diagonal matrix, represented as a vector
131: -  i - INSERT_VALUES or ADD_VALUES

133:    Collective on Mat and Vec

135:    Level: intermediate

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

139: .seealso: MatShift()
140: @*/
141: int MatDiagonalSet(Mat Y,Vec D,InsertMode is)
142: {
143:   int    i,start,end,ierr;

148:   if (Y->ops->diagonalset) {
149:     (*Y->ops->diagonalset)(Y,D,is);
150:   } else {
151:     int    vstart,vend;
152:     PetscScalar *v;
153:     VecGetOwnershipRange(D,&vstart,&vend);
154:     MatGetOwnershipRange(Y,&start,&end);
155:     if (vstart != start || vend != end) {
156:       SETERRQ4(PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %d %d vec %d %d mat",vstart,vend,start,end);
157:     }
158:     VecGetArray(D,&v);
159:     for (i=start; i<end; i++) {
160:       MatSetValues(Y,1,&i,1,&i,v+i-start,is);
161:     }
162:     VecRestoreArray(D,&v);
163:     MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
164:     MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
165:   }
166:   return(0);
167: }

169: #undef __FUNCT__  
171: /*@
172:    MatAYPX - Computes Y = X + a*Y.

174:    Collective on Mat

176:    Input Parameters:
177: +  X,Y - the matrices
178: -  a - the PetscScalar multiplier

180:    Contributed by: Matthew Knepley

182:    Notes:
183:    This routine currently uses the MatAXPY() implementation.

185:    This is slow, if you need it fast send email to petsc-maint@mcs.anl.gov

187:    Level: intermediate

189: .keywords: matrix, add

191: .seealso: MatAXPY()
192:  @*/
193: int MatAYPX(PetscScalar *a,Mat X,Mat Y)
194: {
195:   PetscScalar one = 1.0;
196:   int         mX,mY,nX,nY,ierr;


203:   MatGetSize(X,&mX,&nX);
204:   MatGetSize(X,&mY,&nY);
205:   if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrices: %d %d first %d %d second",mX,mY,nX,nY);

207:   MatScale(a,Y);
208:   MatAXPY(&one,X,Y,DIFFERENT_NONZERO_PATTERN);
209:   return(0);
210: }

212: #undef __FUNCT__  
214: /*@
215:     MatComputeExplicitOperator - Computes the explicit matrix

217:     Collective on Mat

219:     Input Parameter:
220: .   inmat - the matrix

222:     Output Parameter:
223: .   mat - the explict preconditioned operator

225:     Notes:
226:     This computation is done by applying the operators to columns of the 
227:     identity matrix.

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

233:     Level: advanced
234:    
235: .keywords: Mat, compute, explicit, operator

237: @*/
238: int MatComputeExplicitOperator(Mat inmat,Mat *mat)
239: {
240:   Vec      in,out;
241:   int      ierr,i,M,m,size,*rows,start,end;
242:   MPI_Comm comm;
243:   PetscScalar   *array,zero = 0.0,one = 1.0;


249:   comm = inmat->comm;
250:   MPI_Comm_size(comm,&size);

252:   MatGetLocalSize(inmat,&m,0);
253:   MatGetSize(inmat,&M,0);
254:   VecCreateMPI(comm,m,M,&in);
255:   VecDuplicate(in,&out);
256:   VecGetOwnershipRange(in,&start,&end);
257:   PetscMalloc((m+1)*sizeof(int),&rows);
258:   for (i=0; i<m; i++) {rows[i] = start + i;}

260:   if (size == 1) {
261:     MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);
262:   } else {
263:     MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);
264:   }

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

268:     VecSet(&zero,in);
269:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
270:     VecAssemblyBegin(in);
271:     VecAssemblyEnd(in);

273:     MatMult(inmat,in,out);
274: 
275:     VecGetArray(out,&array);
276:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
277:     VecRestoreArray(out,&array);

279:   }
280:   PetscFree(rows);
281:   VecDestroy(out);
282:   VecDestroy(in);
283:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
284:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
285:   return(0);
286: }