Actual source code: axpy.c
petsc-3.4.5 2014-06-29
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(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
157: if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),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
297: .keywords: Mat, compute, explicit, operator
298: @*/
299: PetscErrorCode MatComputeExplicitOperator(Mat inmat,Mat *mat)
300: {
301: Vec in,out;
303: PetscInt i,m,n,M,N,*rows,start,end;
304: MPI_Comm comm;
305: PetscScalar *array,zero = 0.0,one = 1.0;
306: PetscMPIInt size;
312: PetscObjectGetComm((PetscObject)inmat,&comm);
313: MPI_Comm_size(comm,&size);
315: MatGetLocalSize(inmat,&m,&n);
316: MatGetSize(inmat,&M,&N);
317: MatGetVecs(inmat,&in,&out);
318: VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
319: VecGetOwnershipRange(out,&start,&end);
320: PetscMalloc(m*sizeof(PetscInt),&rows);
321: for (i=0; i<m; i++) rows[i] = start + i;
323: MatCreate(comm,mat);
324: MatSetSizes(*mat,m,n,M,N);
325: if (size == 1) {
326: MatSetType(*mat,MATSEQDENSE);
327: MatSeqDenseSetPreallocation(*mat,NULL);
328: } else {
329: MatSetType(*mat,MATMPIAIJ);
330: MatMPIAIJSetPreallocation(*mat,n,NULL,N-n,NULL);
331: }
333: for (i=0; i<N; i++) {
335: VecSet(in,zero);
336: VecSetValues(in,1,&i,&one,INSERT_VALUES);
337: VecAssemblyBegin(in);
338: VecAssemblyEnd(in);
340: MatMult(inmat,in,out);
342: VecGetArray(out,&array);
343: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
344: VecRestoreArray(out,&array);
346: }
347: PetscFree(rows);
348: VecDestroy(&out);
349: VecDestroy(&in);
350: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
351: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
352: return(0);
353: }
355: /* Get the map xtoy which is used by MatAXPY() in the case of SUBSET_NONZERO_PATTERN */
358: PetscErrorCode MatAXPYGetxtoy_Private(PetscInt m,PetscInt *xi,PetscInt *xj,PetscInt *xgarray, PetscInt *yi,PetscInt *yj,PetscInt *ygarray, PetscInt **xtoy)
359: {
361: PetscInt row,i,nz,xcol,ycol,jx,jy,*x2y;
364: PetscMalloc(xi[m]*sizeof(PetscInt),&x2y);
365: i = 0;
366: for (row=0; row<m; row++) {
367: nz = xi[1] - xi[0];
368: jy = 0;
369: for (jx=0; jx<nz; jx++,jy++) {
370: if (xgarray && ygarray) {
371: xcol = xgarray[xj[*xi + jx]];
372: ycol = ygarray[yj[*yi + jy]];
373: } else {
374: xcol = xj[*xi + jx];
375: ycol = yj[*yi + jy]; /* col index for y */
376: }
377: while (ycol < xcol) {
378: jy++;
379: if (ygarray) {
380: ycol = ygarray[yj[*yi + jy]];
381: } else {
382: ycol = yj[*yi + jy];
383: }
384: }
385: if (xcol != ycol) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"X matrix entry (%D,%D) is not in Y matrix",row,ycol);
386: x2y[i++] = *yi + jy;
387: }
388: xi++; yi++;
389: }
390: *xtoy = x2y;
391: return(0);
392: }
396: /*@
397: MatChop - Set all values in the matrix less than the tolerance to zero
399: Input Parameters:
400: + A - The matrix
401: - tol - The zero tolerance
403: Output Parameters:
404: . A - The chopped matrix
406: Level: intermediate
408: .seealso: MatCreate(), MatZeroEntries()
409: @*/
410: PetscErrorCode MatChop(Mat A, PetscReal tol)
411: {
412: PetscScalar *newVals;
413: PetscInt *newCols;
414: PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0;
418: MatGetOwnershipRange(A, &rStart, &rEnd);
419: for (r = rStart; r < rEnd; ++r) {
420: PetscInt ncols;
422: MatGetRow(A, r, &ncols, NULL, NULL);
423: colMax = PetscMax(colMax, ncols);
424: MatRestoreRow(A, r, &ncols, NULL, NULL);
425: }
426: numRows = rEnd - rStart;
427: MPI_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD);
428: PetscMalloc2(colMax,PetscInt,&newCols,colMax,PetscScalar,&newVals);
429: for (r = rStart; r < rStart+maxRows; ++r) {
430: const PetscScalar *vals;
431: const PetscInt *cols;
432: PetscInt ncols, newcols, c;
434: if (r < rEnd) {
435: MatGetRow(A, r, &ncols, &cols, &vals);
436: for (c = 0; c < ncols; ++c) {
437: newCols[c] = cols[c];
438: newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c];
439: }
440: newcols = ncols;
441: MatRestoreRow(A, r, &ncols, &cols, &vals);
442: MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
443: }
444: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
445: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
446: }
447: PetscFree2(newCols,newVals);
448: return(0);
449: }