Actual source code: axpy.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/
6: /*@
7: MatAXPY - Computes Y = a*X + Y.
9: Logically Collective on Mat
11: Input Parameters:
12: + a - the scalar multiplier
13: . X - the first matrix
14: . Y - the second matrix
15: - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN
16: or SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's)
18: Level: intermediate
20: .keywords: matrix, add
22: .seealso: MatAYPX()
23: @*/
24: PetscErrorCode MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)
25: {
27: PetscInt m1,m2,n1,n2;
33: MatGetSize(X,&m1,&n1);
34: MatGetSize(Y,&m2,&n2);
35: if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %D %D %D %D",m1,m2,n1,n2);
37: PetscLogEventBegin(MAT_AXPY,Y,0,0,0);
38: if (Y->ops->axpy) {
39: (*Y->ops->axpy)(Y,a,X,str);
40: } else {
41: MatAXPY_Basic(Y,a,X,str);
42: }
43: PetscLogEventEnd(MAT_AXPY,Y,0,0,0);
44: #if defined(PETSC_HAVE_CUSP)
45: if (Y->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
46: Y->valid_GPU_matrix = PETSC_CUSP_CPU;
47: }
48: #endif
49: return(0);
50: }
54: PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
55: {
56: PetscInt i,start,end,j,ncols,m,n;
57: PetscErrorCode ierr;
58: const PetscInt *row;
59: PetscScalar *val;
60: const PetscScalar *vals;
63: MatGetSize(X,&m,&n);
64: MatGetOwnershipRange(X,&start,&end);
65: if (a == 1.0) {
66: for (i = start; i < end; i++) {
67: MatGetRow(X,i,&ncols,&row,&vals);
68: MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);
69: MatRestoreRow(X,i,&ncols,&row,&vals);
70: }
71: } else {
72: PetscMalloc((n+1)*sizeof(PetscScalar),&val);
73: for (i=start; i<end; i++) {
74: MatGetRow(X,i,&ncols,&row,&vals);
75: for (j=0; j<ncols; j++) {
76: val[j] = a*vals[j];
77: }
78: MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);
79: MatRestoreRow(X,i,&ncols,&row,&vals);
80: }
81: PetscFree(val);
82: }
83: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
84: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
85: return(0);
86: }
90: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str)
91: {
92: PetscInt i,start,end,j,ncols,m,n;
93: PetscErrorCode ierr;
94: const PetscInt *row;
95: PetscScalar *val;
96: const PetscScalar *vals;
99: MatGetSize(X,&m,&n);
100: MatGetOwnershipRange(X,&start,&end);
101: if (a == 1.0) {
102: for (i = start; i < end; i++) {
103: MatGetRow(Y,i,&ncols,&row,&vals);
104: MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
105: MatRestoreRow(Y,i,&ncols,&row,&vals);
107: MatGetRow(X,i,&ncols,&row,&vals);
108: MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
109: MatRestoreRow(X,i,&ncols,&row,&vals);
110: }
111: } else {
112: PetscMalloc((n+1)*sizeof(PetscScalar),&val);
113: for (i=start; i<end; i++) {
114: MatGetRow(Y,i,&ncols,&row,&vals);
115: MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
116: MatRestoreRow(Y,i,&ncols,&row,&vals);
118: MatGetRow(X,i,&ncols,&row,&vals);
119: for (j=0; j<ncols; j++) {
120: val[j] = a*vals[j];
121: }
122: MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);
123: MatRestoreRow(X,i,&ncols,&row,&vals);
124: }
125: PetscFree(val);
126: }
127: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
129: return(0);
130: }
134: /*@
135: MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix.
137: Neighbor-wise Collective on Mat
139: Input Parameters:
140: + Y - the matrices
141: - a - the PetscScalar
143: Level: intermediate
145: .keywords: matrix, add, shift
147: .seealso: MatDiagonalSet()
148: @*/
149: PetscErrorCode MatShift(Mat Y,PetscScalar a)
150: {
152: PetscInt i,start,end;
156: if (!Y->assembled) SETERRQ(((PetscObject)Y)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
157: if (Y->factortype) SETERRQ(((PetscObject)Y)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
158: MatCheckPreallocated(Y,1);
160: if (Y->ops->shift) {
161: (*Y->ops->shift)(Y,a);
162: } else {
163: PetscScalar alpha = a;
164: MatGetOwnershipRange(Y,&start,&end);
165: for (i=start; i<end; i++) {
166: MatSetValues(Y,1,&i,1,&i,&alpha,ADD_VALUES);
167: }
168: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
169: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
170: }
171: #if defined(PETSC_HAVE_CUSP)
172: if (Y->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
173: Y->valid_GPU_matrix = PETSC_CUSP_CPU;
174: }
175: #endif
176: return(0);
177: }
181: PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is)
182: {
184: PetscInt i,start,end,vstart,vend;
185: PetscScalar *v;
188: VecGetOwnershipRange(D,&vstart,&vend);
189: MatGetOwnershipRange(Y,&start,&end);
190: if (vstart != start || vend != end) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %D %D vec %D %D mat",vstart,vend,start,end);
191: VecGetArray(D,&v);
192: for (i=start; i<end; i++) {
193: MatSetValues(Y,1,&i,1,&i,v+i-start,is);
194: }
195: VecRestoreArray(D,&v);
196: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
197: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
198: return(0);
199: }
203: /*@
204: MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
205: that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
206: INSERT_VALUES.
208: Input Parameters:
209: + Y - the input matrix
210: . D - the diagonal matrix, represented as a vector
211: - i - INSERT_VALUES or ADD_VALUES
213: Neighbor-wise Collective on Mat and Vec
215: Level: intermediate
217: .keywords: matrix, add, shift, diagonal
219: .seealso: MatShift()
220: @*/
221: PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is)
222: {
228: if (Y->ops->diagonalset) {
229: (*Y->ops->diagonalset)(Y,D,is);
230: } else {
231: MatDiagonalSet_Default(Y,D,is);
232: }
233: return(0);
234: }
238: /*@
239: MatAYPX - Computes Y = a*Y + X.
241: Logically on Mat
243: Input Parameters:
244: + a - the PetscScalar multiplier
245: . Y - the first matrix
246: . X - the second matrix
247: - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN
249: Level: intermediate
251: .keywords: matrix, add
253: .seealso: MatAXPY()
254: @*/
255: PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str)
256: {
257: PetscScalar one = 1.0;
259: PetscInt mX,mY,nX,nY;
265: MatGetSize(X,&mX,&nX);
266: MatGetSize(X,&mY,&nY);
267: if (mX != mY || nX != nY) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Non conforming matrices: %D %D first %D %D second",mX,mY,nX,nY);
269: MatScale(Y,a);
270: MatAXPY(Y,one,X,str);
271: return(0);
272: }
276: /*@
277: MatComputeExplicitOperator - Computes the explicit matrix
279: Collective on Mat
281: Input Parameter:
282: . inmat - the matrix
284: Output Parameter:
285: . mat - the explict preconditioned operator
287: Notes:
288: This computation is done by applying the operators to columns of the
289: identity matrix.
291: Currently, this routine uses a dense matrix format when 1 processor
292: is used and a sparse format otherwise. This routine is costly in general,
293: and is recommended for use only with relatively small systems.
295: Level: advanced
296:
297: .keywords: Mat, compute, explicit, operator
299: @*/
300: PetscErrorCode MatComputeExplicitOperator(Mat inmat,Mat *mat)
301: {
302: Vec in,out;
304: PetscInt i,m,n,M,N,*rows,start,end;
305: MPI_Comm comm;
306: PetscScalar *array,zero = 0.0,one = 1.0;
307: PetscMPIInt size;
313: comm = ((PetscObject)inmat)->comm;
314: MPI_Comm_size(comm,&size);
316: MatGetLocalSize(inmat,&m,&n);
317: MatGetSize(inmat,&M,&N);
318: MatGetVecs(inmat,&in,&out);
319: VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
320: VecGetOwnershipRange(out,&start,&end);
321: PetscMalloc(m*sizeof(PetscInt),&rows);
322: for (i=0; i<m; i++) {rows[i] = start + i;}
324: MatCreate(comm,mat);
325: MatSetSizes(*mat,m,n,M,N);
326: if (size == 1) {
327: MatSetType(*mat,MATSEQDENSE);
328: MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
329: } else {
330: MatSetType(*mat,MATMPIAIJ);
331: MatMPIAIJSetPreallocation(*mat,n,PETSC_NULL,N-n,PETSC_NULL);
332: }
334: for (i=0; i<N; i++) {
336: VecSet(in,zero);
337: VecSetValues(in,1,&i,&one,INSERT_VALUES);
338: VecAssemblyBegin(in);
339: VecAssemblyEnd(in);
341: MatMult(inmat,in,out);
343: VecGetArray(out,&array);
344: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
345: VecRestoreArray(out,&array);
347: }
348: PetscFree(rows);
349: VecDestroy(&out);
350: VecDestroy(&in);
351: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
352: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
353: return(0);
354: }
356: /* Get the map xtoy which is used by MatAXPY() in the case of SUBSET_NONZERO_PATTERN */
359: PetscErrorCode MatAXPYGetxtoy_Private(PetscInt m,PetscInt *xi,PetscInt *xj,PetscInt *xgarray, PetscInt *yi,PetscInt *yj,PetscInt *ygarray, PetscInt **xtoy)
360: {
362: PetscInt row,i,nz,xcol,ycol,jx,jy,*x2y;
365: PetscMalloc(xi[m]*sizeof(PetscInt),&x2y);
366: i = 0;
367: for (row=0; row<m; row++){
368: nz = xi[1] - xi[0];
369: jy = 0;
370: for (jx=0; jx<nz; jx++,jy++){
371: if (xgarray && ygarray){
372: xcol = xgarray[xj[*xi + jx]];
373: ycol = ygarray[yj[*yi + jy]];
374: } else {
375: xcol = xj[*xi + jx];
376: ycol = yj[*yi + jy]; /* col index for y */
377: }
378: while ( ycol < xcol ) {
379: jy++;
380: if (ygarray){
381: ycol = ygarray[yj[*yi + jy]];
382: } else {
383: ycol = yj[*yi + jy];
384: }
385: }
386: if (xcol != ycol) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"X matrix entry (%D,%D) is not in Y matrix",row,ycol);
387: x2y[i++] = *yi + jy;
388: }
389: xi++; yi++;
390: }
391: *xtoy = x2y;
392: return(0);
393: }