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: }