Actual source code: matutil.c
slepc-3.21.0 2024-03-30
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, G = [ a*A b*B; c*C d*D ].
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: In the case of a parallel matrix, a permuted version of G is returned. The permutation
193: is a perfect shuffle such that the local parts of A, B, C, D remain in the local part of
194: G for the same process.
196: Matrix G must be destroyed by the user.
198: The blocks can be of different type. They can be either ConstantDiagonal, or a standard
199: type such as AIJ, or any other type provided that it supports the MatGetRow operation.
200: The type of the output matrix will be the same as the first block that is not
201: ConstantDiagonal (checked in the A,B,C,D order).
203: Level: developer
205: .seealso: MatCreateNest()
206: @*/
207: PetscErrorCode MatCreateTile(PetscScalar a,Mat A,PetscScalar b,Mat B,PetscScalar c,Mat C,PetscScalar d,Mat D,Mat *G)
208: {
209: PetscInt i,k,M1,M2,N1,N2,M,N,m1,m2,n1,n2,m,n,bs;
210: PetscBool diag[4];
211: Mat block[4] = {A,B,C,D};
212: MatType type[4];
213: PetscMPIInt size;
215: PetscFunctionBegin;
220: PetscCheckSameTypeAndComm(A,2,B,4);
221: PetscCheckSameTypeAndComm(A,2,C,6);
222: PetscCheckSameTypeAndComm(A,2,D,8);
227: PetscAssertPointer(G,9);
229: /* check row 1 */
230: PetscCall(MatGetSize(A,&M1,NULL));
231: PetscCall(MatGetLocalSize(A,&m1,NULL));
232: PetscCall(MatGetSize(B,&M,NULL));
233: PetscCall(MatGetLocalSize(B,&m,NULL));
234: PetscCheck(M==M1 && m==m1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
235: /* check row 2 */
236: PetscCall(MatGetSize(C,&M2,NULL));
237: PetscCall(MatGetLocalSize(C,&m2,NULL));
238: PetscCall(MatGetSize(D,&M,NULL));
239: PetscCall(MatGetLocalSize(D,&m,NULL));
240: PetscCheck(M==M2 && m==m2,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
241: /* check column 1 */
242: PetscCall(MatGetSize(A,NULL,&N1));
243: PetscCall(MatGetLocalSize(A,NULL,&n1));
244: PetscCall(MatGetSize(C,NULL,&N));
245: PetscCall(MatGetLocalSize(C,NULL,&n));
246: PetscCheck(N==N1 && n==n1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
247: /* check column 2 */
248: PetscCall(MatGetSize(B,NULL,&N2));
249: PetscCall(MatGetLocalSize(B,NULL,&n2));
250: PetscCall(MatGetSize(D,NULL,&N));
251: PetscCall(MatGetLocalSize(D,NULL,&n));
252: PetscCheck(N==N2 && n==n2,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible dimensions");
254: /* check matrix types */
255: for (i=0;i<4;i++) {
256: PetscCall(MatGetType(block[i],&type[i]));
257: PetscCall(PetscStrcmp(type[i],MATCONSTANTDIAGONAL,&diag[i]));
258: }
259: for (k=0;k<4;k++) if (!diag[k]) break;
260: PetscCheck(k<4,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not implemented for 4 diagonal blocks");
262: PetscCall(MatGetBlockSize(block[k],&bs));
263: PetscCall(MatCreate(PetscObjectComm((PetscObject)block[k]),G));
264: PetscCall(MatSetSizes(*G,m1+m2,n1+n2,M1+M2,N1+N2));
265: PetscCall(MatSetType(*G,type[k]));
266: PetscCall(MatSetBlockSize(*G,bs));
267: PetscCall(MatSetUp(*G));
269: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)*G),&size));
270: if (size>1) PetscCall(MatCreateTile_MPI(a,A,b,B,c,C,d,D,*G));
271: else PetscCall(MatCreateTile_Seq(a,A,b,B,c,C,d,D,*G));
272: PetscCall(MatAssemblyBegin(*G,MAT_FINAL_ASSEMBLY));
273: PetscCall(MatAssemblyEnd(*G,MAT_FINAL_ASSEMBLY));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*@C
278: MatCreateVecsEmpty - Get vector(s) compatible with the matrix, i.e. with the same
279: parallel layout, but without internal array.
281: Collective
283: Input Parameter:
284: . mat - the matrix
286: Output Parameters:
287: + right - (optional) vector that the matrix can be multiplied against
288: - left - (optional) vector that the matrix vector product can be stored in
290: Note:
291: This is similar to MatCreateVecs(), but the new vectors do not have an internal
292: array, so the intended usage is with VecPlaceArray().
294: Level: developer
296: .seealso: VecDuplicateEmpty()
297: @*/
298: PetscErrorCode MatCreateVecsEmpty(Mat mat,Vec *right,Vec *left)
299: {
300: PetscBool standard,cuda=PETSC_FALSE,skip=PETSC_FALSE;
301: PetscInt M,N,mloc,nloc,rbs,cbs;
302: PetscMPIInt size;
303: Vec v;
305: PetscFunctionBegin;
309: PetscCall(PetscObjectTypeCompareAny((PetscObject)mat,&standard,MATSEQAIJ,MATMPIAIJ,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPISBAIJ,MATSEQDENSE,MATMPIDENSE,""));
310: PetscCall(PetscObjectTypeCompareAny((PetscObject)mat,&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,""));
311: if (!standard && !cuda) {
312: PetscCall(MatCreateVecs(mat,right,left));
313: v = right? *right: *left;
314: if (v) {
315: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,""));
316: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,""));
317: }
318: if (!standard && !cuda) skip = PETSC_TRUE;
319: else {
320: if (right) PetscCall(VecDestroy(right));
321: if (left) PetscCall(VecDestroy(left));
322: }
323: }
324: if (!skip) {
325: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size));
326: PetscCall(MatGetLocalSize(mat,&mloc,&nloc));
327: PetscCall(MatGetSize(mat,&M,&N));
328: PetscCall(MatGetBlockSizes(mat,&rbs,&cbs));
329: if (right) {
330: if (cuda) {
331: #if defined(PETSC_HAVE_CUDA)
332: if (size>1) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)mat),cbs,nloc,N,NULL,right));
333: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)mat),cbs,N,NULL,right));
334: #endif
335: } else {
336: if (size>1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),cbs,nloc,N,NULL,right));
337: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)mat),cbs,N,NULL,right));
338: }
339: }
340: if (left) {
341: if (cuda) {
342: #if defined(PETSC_HAVE_CUDA)
343: if (size>1) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)mat),rbs,mloc,M,NULL,left));
344: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)mat),rbs,M,NULL,left));
345: #endif
346: } else {
347: if (size>1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),rbs,mloc,M,NULL,left));
348: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)mat),rbs,M,NULL,left));
349: }
350: }
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@C
356: MatNormEstimate - Estimate the 2-norm of a matrix.
358: Collective
360: Input Parameters:
361: + A - the matrix
362: . vrn - random vector with normally distributed entries (can be NULL)
363: - w - workspace vector (can be NULL)
365: Output Parameter:
366: . nrm - the norm estimate
368: Notes:
369: Does not need access to the matrix entries, just performs a matrix-vector product.
370: Based on work by I. Ipsen and coworkers https://ipsen.math.ncsu.edu/ps/slides_ima.pdf
372: The input vector vrn must have unit 2-norm.
373: If vrn is NULL, then it is created internally and filled with VecSetRandomNormal().
375: Level: developer
377: .seealso: VecSetRandomNormal()
378: @*/
379: PetscErrorCode MatNormEstimate(Mat A,Vec vrn,Vec w,PetscReal *nrm)
380: {
381: PetscInt n;
382: Vec vv=NULL,ww=NULL;
384: PetscFunctionBegin;
389: PetscAssertPointer(nrm,4);
391: if (!vrn) {
392: PetscCall(MatCreateVecs(A,&vv,NULL));
393: vrn = vv;
394: PetscCall(VecSetRandomNormal(vv,NULL,NULL,NULL));
395: PetscCall(VecNormalize(vv,NULL));
396: }
397: if (!w) {
398: PetscCall(MatCreateVecs(A,&ww,NULL));
399: w = ww;
400: }
402: PetscCall(MatGetSize(A,&n,NULL));
403: PetscCall(MatMult(A,vrn,w));
404: PetscCall(VecNorm(w,NORM_2,nrm));
405: *nrm *= PetscSqrtReal((PetscReal)n);
407: PetscCall(VecDestroy(&vv));
408: PetscCall(VecDestroy(&ww));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }