Actual source code: axpy.c
petsc-3.6.1 2015-08-06
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: #if defined(PETSC_HAVE_VIENNACL)
50: if (Y->valid_GPU_matrix != PETSC_VIENNACL_UNALLOCATED) {
51: Y->valid_GPU_matrix = PETSC_VIENNACL_CPU;
52: }
53: #endif
54: return(0);
55: }
59: PetscErrorCode MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)
60: {
61: PetscInt i,start,end,j,ncols,m,n;
62: PetscErrorCode ierr;
63: const PetscInt *row;
64: PetscScalar *val;
65: const PetscScalar *vals;
68: MatGetSize(X,&m,&n);
69: MatGetOwnershipRange(X,&start,&end);
70: if (a == 1.0) {
71: for (i = start; i < end; i++) {
72: MatGetRow(X,i,&ncols,&row,&vals);
73: MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);
74: MatRestoreRow(X,i,&ncols,&row,&vals);
75: }
76: } else {
77: PetscMalloc1(n+1,&val);
78: for (i=start; i<end; i++) {
79: MatGetRow(X,i,&ncols,&row,&vals);
80: for (j=0; j<ncols; j++) {
81: val[j] = a*vals[j];
82: }
83: MatSetValues(Y,1,&i,ncols,row,val,ADD_VALUES);
84: MatRestoreRow(X,i,&ncols,&row,&vals);
85: }
86: PetscFree(val);
87: }
88: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
90: return(0);
91: }
95: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str)
96: {
97: PetscInt i,start,end,j,ncols,m,n;
98: PetscErrorCode ierr;
99: const PetscInt *row;
100: PetscScalar *val;
101: const PetscScalar *vals;
104: MatGetSize(X,&m,&n);
105: MatGetOwnershipRange(X,&start,&end);
106: if (a == 1.0) {
107: for (i = start; i < end; i++) {
108: MatGetRow(Y,i,&ncols,&row,&vals);
109: MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
110: MatRestoreRow(Y,i,&ncols,&row,&vals);
112: MatGetRow(X,i,&ncols,&row,&vals);
113: MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
114: MatRestoreRow(X,i,&ncols,&row,&vals);
115: }
116: } else {
117: PetscMalloc1(n+1,&val);
118: for (i=start; i<end; i++) {
119: MatGetRow(Y,i,&ncols,&row,&vals);
120: MatSetValues(B,1,&i,ncols,row,vals,ADD_VALUES);
121: MatRestoreRow(Y,i,&ncols,&row,&vals);
123: MatGetRow(X,i,&ncols,&row,&vals);
124: for (j=0; j<ncols; j++) {
125: val[j] = a*vals[j];
126: }
127: MatSetValues(B,1,&i,ncols,row,val,ADD_VALUES);
128: MatRestoreRow(X,i,&ncols,&row,&vals);
129: }
130: PetscFree(val);
131: }
132: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
133: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
134: return(0);
135: }
139: /*@
140: MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix.
142: Neighbor-wise Collective on Mat
144: Input Parameters:
145: + Y - the matrices
146: - a - the PetscScalar
148: Level: intermediate
150: .keywords: matrix, add, shift
152: .seealso: MatDiagonalSet()
153: @*/
154: PetscErrorCode MatShift(Mat Y,PetscScalar a)
155: {
160: if (!Y->assembled) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
161: if (Y->factortype) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
162: MatCheckPreallocated(Y,1);
164: (*Y->ops->shift)(Y,a);
166: #if defined(PETSC_HAVE_CUSP)
167: if (Y->valid_GPU_matrix != PETSC_CUSP_UNALLOCATED) {
168: Y->valid_GPU_matrix = PETSC_CUSP_CPU;
169: }
170: #endif
171: #if defined(PETSC_HAVE_VIENNACL)
172: if (Y->valid_GPU_matrix != PETSC_VIENNACL_UNALLOCATED) {
173: Y->valid_GPU_matrix = PETSC_VIENNACL_CPU;
174: }
175: #endif
176: return(0);
177: }
181: PetscErrorCode MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is)
182: {
184: PetscInt i,start,end;
185: PetscScalar *v;
188: MatGetOwnershipRange(Y,&start,&end);
189: VecGetArray(D,&v);
190: for (i=start; i<end; i++) {
191: MatSetValues(Y,1,&i,1,&i,v+i-start,is);
192: }
193: VecRestoreArray(D,&v);
194: MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
195: MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
196: return(0);
197: }
201: /*@
202: MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix
203: that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is
204: INSERT_VALUES.
206: Input Parameters:
207: + Y - the input matrix
208: . D - the diagonal matrix, represented as a vector
209: - i - INSERT_VALUES or ADD_VALUES
211: Neighbor-wise Collective on Mat and Vec
213: Level: intermediate
215: .keywords: matrix, add, shift, diagonal
217: .seealso: MatShift()
218: @*/
219: PetscErrorCode MatDiagonalSet(Mat Y,Vec D,InsertMode is)
220: {
222: PetscInt matlocal,veclocal;
227: MatGetLocalSize(Y,&matlocal,NULL);
228: VecGetLocalSize(D,&veclocal);
229: if (matlocal != veclocal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number local rows of matrix %D does not match that of vector for diagonal %D",matlocal,veclocal);
230: if (Y->ops->diagonalset) {
231: (*Y->ops->diagonalset)(Y,D,is);
232: } else {
233: MatDiagonalSet_Default(Y,D,is);
234: }
235: return(0);
236: }
240: /*@
241: MatAYPX - Computes Y = a*Y + X.
243: Logically on Mat
245: Input Parameters:
246: + a - the PetscScalar multiplier
247: . Y - the first matrix
248: . X - the second matrix
249: - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN
251: Level: intermediate
253: .keywords: matrix, add
255: .seealso: MatAXPY()
256: @*/
257: PetscErrorCode MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str)
258: {
259: PetscScalar one = 1.0;
261: PetscInt mX,mY,nX,nY;
267: MatGetSize(X,&mX,&nX);
268: MatGetSize(X,&mY,&nY);
269: 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);
271: MatScale(Y,a);
272: MatAXPY(Y,one,X,str);
273: return(0);
274: }
278: /*@
279: MatComputeExplicitOperator - Computes the explicit matrix
281: Collective on Mat
283: Input Parameter:
284: . inmat - the matrix
286: Output Parameter:
287: . mat - the explict preconditioned operator
289: Notes:
290: This computation is done by applying the operators to columns of the
291: identity matrix.
293: Currently, this routine uses a dense matrix format when 1 processor
294: is used and a sparse format otherwise. This routine is costly in general,
295: and is recommended for use only with relatively small systems.
297: Level: advanced
299: .keywords: Mat, compute, explicit, operator
300: @*/
301: PetscErrorCode MatComputeExplicitOperator(Mat inmat,Mat *mat)
302: {
303: Vec in,out;
305: PetscInt i,m,n,M,N,*rows,start,end;
306: MPI_Comm comm;
307: PetscScalar *array,zero = 0.0,one = 1.0;
308: PetscMPIInt size;
314: PetscObjectGetComm((PetscObject)inmat,&comm);
315: MPI_Comm_size(comm,&size);
317: MatGetLocalSize(inmat,&m,&n);
318: MatGetSize(inmat,&M,&N);
319: MatCreateVecs(inmat,&in,&out);
320: VecSetOption(in,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
321: VecGetOwnershipRange(out,&start,&end);
322: PetscMalloc1(m,&rows);
323: for (i=0; i<m; i++) rows[i] = start + i;
325: MatCreate(comm,mat);
326: MatSetSizes(*mat,m,n,M,N);
327: if (size == 1) {
328: MatSetType(*mat,MATSEQDENSE);
329: MatSeqDenseSetPreallocation(*mat,NULL);
330: } else {
331: MatSetType(*mat,MATMPIAIJ);
332: MatMPIAIJSetPreallocation(*mat,n,NULL,N-n,NULL);
333: }
335: for (i=0; i<N; i++) {
337: VecSet(in,zero);
338: VecSetValues(in,1,&i,&one,INSERT_VALUES);
339: VecAssemblyBegin(in);
340: VecAssemblyEnd(in);
342: MatMult(inmat,in,out);
344: VecGetArray(out,&array);
345: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
346: VecRestoreArray(out,&array);
348: }
349: PetscFree(rows);
350: VecDestroy(&out);
351: VecDestroy(&in);
352: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
353: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
354: return(0);
355: }
359: /*@
360: MatChop - Set all values in the matrix less than the tolerance to zero
362: Input Parameters:
363: + A - The matrix
364: - tol - The zero tolerance
366: Output Parameters:
367: . A - The chopped matrix
369: Level: intermediate
371: .seealso: MatCreate(), MatZeroEntries()
372: @*/
373: PetscErrorCode MatChop(Mat A, PetscReal tol)
374: {
375: PetscScalar *newVals;
376: PetscInt *newCols;
377: PetscInt rStart, rEnd, numRows, maxRows, r, colMax = 0;
381: MatGetOwnershipRange(A, &rStart, &rEnd);
382: for (r = rStart; r < rEnd; ++r) {
383: PetscInt ncols;
385: MatGetRow(A, r, &ncols, NULL, NULL);
386: colMax = PetscMax(colMax, ncols);
387: MatRestoreRow(A, r, &ncols, NULL, NULL);
388: }
389: numRows = rEnd - rStart;
390: MPI_Allreduce(&numRows, &maxRows, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)A));
391: PetscMalloc2(colMax,&newCols,colMax,&newVals);
392: for (r = rStart; r < rStart+maxRows; ++r) {
393: const PetscScalar *vals;
394: const PetscInt *cols;
395: PetscInt ncols, newcols, c;
397: if (r < rEnd) {
398: MatGetRow(A, r, &ncols, &cols, &vals);
399: for (c = 0; c < ncols; ++c) {
400: newCols[c] = cols[c];
401: newVals[c] = PetscAbsScalar(vals[c]) < tol ? 0.0 : vals[c];
402: }
403: newcols = ncols;
404: MatRestoreRow(A, r, &ncols, &cols, &vals);
405: MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
406: }
407: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
408: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
409: }
410: PetscFree2(newCols,newVals);
411: return(0);
412: }