Actual source code: matutil.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/slepcimpl.h>
13: static PetscErrorCode MatCreateTile_Seq(PetscScalar a,Mat A,PetscScalar b,Mat B,PetscScalar c,Mat C,PetscScalar d,Mat D,Mat G)
14: {
15: PetscInt i,j,M1,M2,N1,N2,ncols,*scols;
16: PetscScalar *svals,*buf;
17: const PetscInt *cols;
18: const PetscScalar *vals;
20: PetscFunctionBegin;
21: PetscCall(MatGetSize(A,&M1,&N1));
22: PetscCall(MatGetSize(D,&M2,&N2));
24: PetscCall(PetscMalloc2(PetscMax(N1,N2),&buf,PetscMax(N1,N2),&scols));
25: /* Transfer A */
26: if (a!=0.0) {
27: for (i=0;i<M1;i++) {
28: PetscCall(MatGetRow(A,i,&ncols,&cols,&vals));
29: if (a!=1.0) {
30: svals=buf;
31: for (j=0;j<ncols;j++) svals[j] = vals[j]*a;
32: } else svals=(PetscScalar*)vals;
33: PetscCall(MatSetValues(G,1,&i,ncols,cols,svals,INSERT_VALUES));
34: PetscCall(MatRestoreRow(A,i,&ncols,&cols,&vals));
35: }
36: }
37: /* Transfer B */
38: if (b!=0.0) {
39: for (i=0;i<M1;i++) {
40: PetscCall(MatGetRow(B,i,&ncols,&cols,&vals));
41: if (b!=1.0) {
42: svals=buf;
43: for (j=0;j<ncols;j++) svals[j] = vals[j]*b;
44: } else svals=(PetscScalar*)vals;
45: for (j=0;j<ncols;j++) scols[j] = cols[j]+N1;
46: PetscCall(MatSetValues(G,1,&i,ncols,scols,svals,INSERT_VALUES));
47: PetscCall(MatRestoreRow(B,i,&ncols,&cols,&vals));
48: }
49: }
50: /* Transfer C */
51: if (c!=0.0) {
52: for (i=0;i<M2;i++) {
53: PetscCall(MatGetRow(C,i,&ncols,&cols,&vals));
54: if (c!=1.0) {
55: svals=buf;
56: for (j=0;j<ncols;j++) svals[j] = vals[j]*c;
57: } else svals=(PetscScalar*)vals;
58: j = i+M1;
59: PetscCall(MatSetValues(G,1,&j,ncols,cols,svals,INSERT_VALUES));
60: PetscCall(MatRestoreRow(C,i,&ncols,&cols,&vals));
61: }
62: }
63: /* Transfer D */
64: if (d!=0.0) {
65: for (i=0;i<M2;i++) {
66: PetscCall(MatGetRow(D,i,&ncols,&cols,&vals));
67: if (d!=1.0) {
68: svals=buf;
69: for (j=0;j<ncols;j++) svals[j] = vals[j]*d;
70: } else svals=(PetscScalar*)vals;
71: for (j=0;j<ncols;j++) scols[j] = cols[j]+N1;
72: j = i+M1;
73: PetscCall(MatSetValues(G,1,&j,ncols,scols,svals,INSERT_VALUES));
74: PetscCall(MatRestoreRow(D,i,&ncols,&cols,&vals));
75: }
76: }
77: PetscCall(PetscFree2(buf,scols));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: static PetscErrorCode MatCreateTile_MPI(PetscScalar a,Mat A,PetscScalar b,Mat B,PetscScalar c,Mat C,PetscScalar d,Mat D,Mat G)
82: {
83: PetscMPIInt np;
84: PetscInt p,i,j,N1,N2,m1,m2,*map1,*map2;
85: PetscInt ncols,*scols,start,gstart;
86: PetscScalar *svals,*buf;
87: const PetscInt *cols,*mapptr1,*mapptr2;
88: const PetscScalar *vals;
90: PetscFunctionBegin;
91: PetscCall(MatGetSize(A,NULL,&N1));
92: PetscCall(MatGetLocalSize(A,&m1,NULL));
93: PetscCall(MatGetSize(D,NULL,&N2));
94: PetscCall(MatGetLocalSize(D,&m2,NULL));
96: /* Create mappings */
97: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)G),&np));
98: PetscCall(MatGetOwnershipRangesColumn(A,&mapptr1));
99: PetscCall(MatGetOwnershipRangesColumn(B,&mapptr2));
100: PetscCall(PetscMalloc4(PetscMax(N1,N2),&buf,PetscMax(N1,N2),&scols,N1,&map1,N2,&map2));
101: for (p=0;p<np;p++) {
102: for (i=mapptr1[p];i<mapptr1[p+1];i++) map1[i] = i+mapptr2[p];
103: }
104: for (p=0;p<np;p++) {
105: for (i=mapptr2[p];i<mapptr2[p+1];i++) map2[i] = i+mapptr1[p+1];
106: }
107: PetscCall(MatGetOwnershipRange(G,&gstart,NULL));
109: /* Transfer A */
110: if (a!=0.0) {
111: PetscCall(MatGetOwnershipRange(A,&start,NULL));
112: for (i=0;i<m1;i++) {
113: PetscCall(MatGetRow(A,i+start,&ncols,&cols,&vals));
114: if (a!=1.0) {
115: svals=buf;
116: for (j=0;j<ncols;j++) svals[j] = vals[j]*a;
117: } else svals=(PetscScalar*)vals;
118: for (j=0;j<ncols;j++) scols[j] = map1[cols[j]];
119: j = gstart+i;
120: PetscCall(MatSetValues(G,1,&j,ncols,scols,svals,INSERT_VALUES));
121: PetscCall(MatRestoreRow(A,i+start,&ncols,&cols,&vals));
122: }
123: }
124: /* Transfer B */
125: if (b!=0.0) {
126: PetscCall(MatGetOwnershipRange(B,&start,NULL));
127: for (i=0;i<m1;i++) {
128: PetscCall(MatGetRow(B,i+start,&ncols,&cols,&vals));
129: if (b!=1.0) {
130: svals=buf;
131: for (j=0;j<ncols;j++) svals[j] = vals[j]*b;
132: } else svals=(PetscScalar*)vals;
133: for (j=0;j<ncols;j++) scols[j] = map2[cols[j]];
134: j = gstart+i;
135: PetscCall(MatSetValues(G,1,&j,ncols,scols,svals,INSERT_VALUES));
136: PetscCall(MatRestoreRow(B,i+start,&ncols,&cols,&vals));
137: }
138: }
139: /* Transfer C */
140: if (c!=0.0) {
141: PetscCall(MatGetOwnershipRange(C,&start,NULL));
142: for (i=0;i<m2;i++) {
143: PetscCall(MatGetRow(C,i+start,&ncols,&cols,&vals));
144: if (c!=1.0) {
145: svals=buf;
146: for (j=0;j<ncols;j++) svals[j] = vals[j]*c;
147: } else svals=(PetscScalar*)vals;
148: for (j=0;j<ncols;j++) scols[j] = map1[cols[j]];
149: j = gstart+m1+i;
150: PetscCall(MatSetValues(G,1,&j,ncols,scols,svals,INSERT_VALUES));
151: PetscCall(MatRestoreRow(C,i+start,&ncols,&cols,&vals));
152: }
153: }
154: /* Transfer D */
155: if (d!=0.0) {
156: PetscCall(MatGetOwnershipRange(D,&start,NULL));
157: for (i=0;i<m2;i++) {
158: PetscCall(MatGetRow(D,i+start,&ncols,&cols,&vals));
159: if (d!=1.0) {
160: svals=buf;
161: for (j=0;j<ncols;j++) svals[j] = vals[j]*d;
162: } else svals=(PetscScalar*)vals;
163: for (j=0;j<ncols;j++) scols[j] = map2[cols[j]];
164: j = gstart+m1+i;
165: PetscCall(MatSetValues(G,1,&j,ncols,scols,svals,INSERT_VALUES));
166: PetscCall(MatRestoreRow(D,i+start,&ncols,&cols,&vals));
167: }
168: }
169: PetscCall(PetscFree4(buf,scols,map1,map2));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@
174: MatCreateTile - Explicitly build a matrix from four blocks.
176: Collective
178: Input Parameters:
179: + A - matrix for top-left block
180: . a - scaling factor for block A
181: . B - matrix for top-right block
182: . b - scaling factor for block B
183: . C - matrix for bottom-left block
184: . c - scaling factor for block C
185: . D - matrix for bottom-right block
186: - d - scaling factor for block D
188: Output Parameter:
189: . G - the resulting matrix
191: Notes:
192: The resulting matrix is built as
193: $$G = \begin{bmatrix} aA & bB\\ cC & dD \end{bmatrix}.$$
195: In the case of a parallel matrix, a permuted version of `G` is returned. The permutation
196: is a perfect shuffle such that the local parts of `A`, `B`, `C`, `D` remain in the local part of
197: `G` for the same process.
199: Matrix `G` must be destroyed by the user.
201: The blocks can be of different type. They can be either `MATCONSTANTDIAGONAL`, or a standard
202: type such as `MATAIJ`, or any other type provided that it supports the `MatGetRow()` operation.
203: The type of the output matrix will be the same as the first block that is not
204: `MATCONSTANTDIAGONAL` (checked in the `A`,`B`,`C`,`D` order).
206: Level: developer
208: .seealso: `MatCreateNest()`
209: @*/
210: PetscErrorCode MatCreateTile(PetscScalar a,Mat A,PetscScalar b,Mat B,PetscScalar c,Mat C,PetscScalar d,Mat D,Mat *G)
211: {
212: PetscInt i,k,M1,M2,N1,N2,M,N,m1,m2,n1,n2,m,n,bs;
213: PetscBool diag[4];
214: Mat block[4] = {A,B,C,D};
215: MatType type[4];
216: PetscMPIInt size;
218: PetscFunctionBegin;
223: PetscCheckSameTypeAndComm(A,2,B,4);
224: PetscCheckSameTypeAndComm(A,2,C,6);
225: PetscCheckSameTypeAndComm(A,2,D,8);
230: PetscAssertPointer(G,9);
232: /* check row 1 */
233: PetscCall(MatGetSize(A,&M1,NULL));
234: PetscCall(MatGetLocalSize(A,&m1,NULL));
235: PetscCall(MatGetSize(B,&M,NULL));
236: PetscCall(MatGetLocalSize(B,&m,NULL));
237: PetscCheck(M==M1 && m==m1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
238: /* check row 2 */
239: PetscCall(MatGetSize(C,&M2,NULL));
240: PetscCall(MatGetLocalSize(C,&m2,NULL));
241: PetscCall(MatGetSize(D,&M,NULL));
242: PetscCall(MatGetLocalSize(D,&m,NULL));
243: PetscCheck(M==M2 && m==m2,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
244: /* check column 1 */
245: PetscCall(MatGetSize(A,NULL,&N1));
246: PetscCall(MatGetLocalSize(A,NULL,&n1));
247: PetscCall(MatGetSize(C,NULL,&N));
248: PetscCall(MatGetLocalSize(C,NULL,&n));
249: PetscCheck(N==N1 && n==n1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
250: /* check column 2 */
251: PetscCall(MatGetSize(B,NULL,&N2));
252: PetscCall(MatGetLocalSize(B,NULL,&n2));
253: PetscCall(MatGetSize(D,NULL,&N));
254: PetscCall(MatGetLocalSize(D,NULL,&n));
255: PetscCheck(N==N2 && n==n2,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
257: /* check matrix types */
258: for (i=0;i<4;i++) {
259: PetscCall(MatGetType(block[i],&type[i]));
260: PetscCall(PetscStrcmp(type[i],MATCONSTANTDIAGONAL,&diag[i]));
261: }
262: for (k=0;k<4;k++) if (!diag[k]) break;
263: PetscCheck(k<4,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented for 4 diagonal blocks");
265: PetscCall(MatGetBlockSize(block[k],&bs));
266: PetscCall(MatCreate(PetscObjectComm((PetscObject)block[k]),G));
267: PetscCall(MatSetSizes(*G,m1+m2,n1+n2,M1+M2,N1+N2));
268: PetscCall(MatSetType(*G,type[k]));
269: PetscCall(MatSetBlockSize(*G,bs));
270: PetscCall(MatSetUp(*G));
272: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)*G),&size));
273: if (size>1) PetscCall(MatCreateTile_MPI(a,A,b,B,c,C,d,D,*G));
274: else PetscCall(MatCreateTile_Seq(a,A,b,B,c,C,d,D,*G));
275: PetscCall(MatAssemblyBegin(*G,MAT_FINAL_ASSEMBLY));
276: PetscCall(MatAssemblyEnd(*G,MAT_FINAL_ASSEMBLY));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@
281: MatCreateVecsEmpty - Get vector(s) compatible with the matrix, i.e., with the same
282: parallel layout, but without internal array.
284: Collective
286: Input Parameter:
287: . mat - the matrix
289: Output Parameters:
290: + right - (optional) vector that the matrix can be multiplied against
291: - left - (optional) vector that the matrix vector product can be stored in
293: Note:
294: This is similar to `MatCreateVecs()`, but the new vectors do not have an internal
295: array, so the intended usage is with `VecPlaceArray()`.
297: Level: developer
299: .seealso: `MatCreateVecs()`, `VecDuplicateEmpty()`
300: @*/
301: PetscErrorCode MatCreateVecsEmpty(Mat mat,Vec *right,Vec *left)
302: {
303: PetscBool standard,cuda=PETSC_FALSE,hip=PETSC_FALSE,skip=PETSC_FALSE;
304: PetscInt M,N,mloc,nloc,rbs,cbs;
305: PetscMPIInt size;
306: Vec v;
308: PetscFunctionBegin;
312: PetscCall(PetscObjectTypeCompareAny((PetscObject)mat,&standard,MATSEQAIJ,MATMPIAIJ,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPISBAIJ,MATSEQDENSE,MATMPIDENSE,""));
313: PetscCall(PetscObjectTypeCompareAny((PetscObject)mat,&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
314: PetscCall(PetscObjectTypeCompareAny((PetscObject)mat,&hip,MATSEQAIJHIPSPARSE,MATMPIAIJHIPSPARSE,""));
315: if (!standard && !cuda && !hip) {
316: PetscCall(MatCreateVecs(mat,right,left));
317: v = right? *right: *left;
318: if (v) {
319: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,""));
320: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,""));
321: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&hip,VECSEQHIP,VECMPIHIP,""));
322: }
323: if (!standard && !cuda && !hip) skip = PETSC_TRUE;
324: else {
325: if (right) PetscCall(VecDestroy(right));
326: if (left) PetscCall(VecDestroy(left));
327: }
328: }
329: if (!skip) {
330: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size));
331: PetscCall(MatGetLocalSize(mat,&mloc,&nloc));
332: PetscCall(MatGetSize(mat,&M,&N));
333: PetscCall(MatGetBlockSizes(mat,&rbs,&cbs));
334: if (right) {
335: if (cuda) {
336: #if defined(PETSC_HAVE_CUDA)
337: if (size>1) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)mat),cbs,nloc,N,NULL,right));
338: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)mat),cbs,N,NULL,right));
339: #endif
340: } else if (hip) {
341: #if defined(PETSC_HAVE_HIP)
342: if (size>1) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)mat),cbs,nloc,N,NULL,right));
343: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)mat),cbs,N,NULL,right));
344: #endif
345: } else {
346: if (size>1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),cbs,nloc,N,NULL,right));
347: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)mat),cbs,N,NULL,right));
348: }
349: }
350: if (left) {
351: if (cuda) {
352: #if defined(PETSC_HAVE_CUDA)
353: if (size>1) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)mat),rbs,mloc,M,NULL,left));
354: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)mat),rbs,M,NULL,left));
355: #endif
356: } else if (hip) {
357: #if defined(PETSC_HAVE_HIP)
358: if (size>1) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)mat),rbs,mloc,M,NULL,left));
359: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)mat),rbs,M,NULL,left));
360: #endif
361: } else {
362: if (size>1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),rbs,mloc,M,NULL,left));
363: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)mat),rbs,M,NULL,left));
364: }
365: }
366: }
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: /*@
371: MatNormEstimate - Estimate the 2-norm of a matrix.
373: Collective
375: Input Parameters:
376: + A - the matrix
377: . vrn - random vector with normally distributed entries (can be `NULL`)
378: - w - workspace vector (can be `NULL`)
380: Output Parameter:
381: . nrm - the norm estimate
383: Notes:
384: Does not need access to the matrix entries, just performs a matrix-vector product.
385: Based on work by I. Ipsen and coworkers <https://ipsen.math.ncsu.edu/ps/slides_ima.pdf>.
387: The input vector `vrn` must have unit 2-norm.
388: If `vrn` is `NULL`, then it is created internally and filled with `VecSetRandomNormal()`.
390: Level: developer
392: .seealso: `VecSetRandomNormal()`
393: @*/
394: PetscErrorCode MatNormEstimate(Mat A,Vec vrn,Vec w,PetscReal *nrm)
395: {
396: PetscInt n;
397: Vec vv=NULL,ww=NULL;
399: PetscFunctionBegin;
404: PetscAssertPointer(nrm,4);
406: if (!vrn) {
407: PetscCall(MatCreateVecs(A,&vv,NULL));
408: vrn = vv;
409: PetscCall(VecSetRandomNormal(vv,NULL,NULL,NULL));
410: PetscCall(VecNormalize(vv,NULL));
411: }
412: if (!w) {
413: PetscCall(MatCreateVecs(A,&ww,NULL));
414: w = ww;
415: }
417: PetscCall(MatGetSize(A,&n,NULL));
418: PetscCall(MatMult(A,vrn,w));
419: PetscCall(VecNorm(w,NORM_2,nrm));
420: *nrm *= PetscSqrtReal((PetscReal)n);
422: PetscCall(VecDestroy(&vv));
423: PetscCall(VecDestroy(&ww));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }