Actual source code: ex192.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Tests MatSolve() and MatMatSolve() with MUMPS or MKL_PARDISO sequential solvers in Schur complement mode.\n\
3: Example: mpiexec -n 1 ./ex192 -f <matrix binary file> -nrhs 4 -symmetric_solve -hermitian_solve -schur_ratio 0.3\n\n";
5: #include <petscmat.h>
7: int main(int argc,char **args)
8: {
9: Mat A,RHS,C,F,X,S;
10: Vec u,x,b;
11: Vec xschur,bschur,uschur;
12: IS is_schur;
14: PetscMPIInt size;
15: PetscInt isolver=0,size_schur,m,n,nfact,nsolve,nrhs;
16: PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON;
17: PetscRandom rand;
18: PetscBool data_provided,herm,symm,use_lu;
19: PetscReal sratio = 5.1/12.;
20: PetscViewer fd; /* viewer */
21: char solver[256];
22: char file[PETSC_MAX_PATH_LEN]; /* input file name */
24: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
25: MPI_Comm_size(PETSC_COMM_WORLD, &size);
26: if (size > 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor test");
27: /* Determine which type of solver we want to test for */
28: herm = PETSC_FALSE;
29: symm = PETSC_FALSE;
30: PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);
31: PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);
32: if (herm) symm = PETSC_TRUE;
34: /* Determine file from which we read the matrix A */
35: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&data_provided);
36: if (!data_provided) { /* get matrices from PETSc distribution */
37: PetscStrncpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/",sizeof(file));
38: if (symm) {
39: #if defined (PETSC_USE_COMPLEX)
40: PetscStrlcat(file,"hpd-complex-",sizeof(file));
41: #else
42: PetscStrlcat(file,"spd-real-",sizeof(file));
43: #endif
44: } else {
45: #if defined (PETSC_USE_COMPLEX)
46: PetscStrlcat(file,"nh-complex-",sizeof(file));
47: #else
48: PetscStrlcat(file,"ns-real-",sizeof(file));
49: #endif
50: }
51: #if defined(PETSC_USE_64BIT_INDICES)
52: PetscStrlcat(file,"int64-",sizeof(file));
53: #else
54: PetscStrlcat(file,"int32-",sizeof(file));
55: #endif
56: #if defined (PETSC_USE_REAL_SINGLE)
57: PetscStrlcat(file,"float32",sizeof(file));
58: #else
59: PetscStrlcat(file,"float64",sizeof(file));
60: #endif
61: }
62: /* Load matrix A */
63: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
64: MatCreate(PETSC_COMM_WORLD,&A);
65: MatLoad(A,fd);
66: PetscViewerDestroy(&fd);
67: MatGetSize(A,&m,&n);
68: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
70: /* Create dense matrix C and X; C holds true solution with identical colums */
71: nrhs = 2;
72: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
73: MatCreate(PETSC_COMM_WORLD,&C);
74: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
75: MatSetType(C,MATDENSE);
76: MatSetFromOptions(C);
77: MatSetUp(C);
79: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
80: PetscRandomSetFromOptions(rand);
81: MatSetRandom(C,rand);
82: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
84: /* Create vectors */
85: VecCreate(PETSC_COMM_WORLD,&x);
86: VecSetSizes(x,n,PETSC_DECIDE);
87: VecSetFromOptions(x);
88: VecDuplicate(x,&b);
89: VecDuplicate(x,&u); /* save the true solution */
91: PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);
92: switch (isolver) {
93: #if defined(PETSC_HAVE_MUMPS)
94: case 0:
95: PetscStrcpy(solver,MATSOLVERMUMPS);
96: break;
97: #endif
98: #if defined(PETSC_HAVE_MKL_PARDISO)
99: case 1:
100: PetscStrcpy(solver,MATSOLVERMKL_PARDISO);
101: break;
102: #endif
103: default:
104: PetscStrcpy(solver,MATSOLVERPETSC);
105: break;
106: }
108: #if defined (PETSC_USE_COMPLEX)
109: if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
110: PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
111: PetscScalar val = -1.0;
112: val = val + im;
113: MatSetValue(A,1,0,val,INSERT_VALUES);
114: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
116: }
117: #endif
119: PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);
120: if (sratio < 0. || sratio > 1.) {
121: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
122: }
123: size_schur = (PetscInt)(sratio*m);
125: PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %D, sym %d, herm %d, size schur %D, size mat %D\n",solver,nrhs,symm,herm,size_schur,m);
127: /* Test LU/Cholesky Factorization */
128: use_lu = PETSC_FALSE;
129: if (!symm) use_lu = PETSC_TRUE;
130: #if defined (PETSC_USE_COMPLEX)
131: if (isolver == 1) use_lu = PETSC_TRUE;
132: #endif
134: if (herm && !use_lu) { /* test also conversion routines inside the solver packages */
135: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
136: MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);
137: }
140: if (use_lu) {
141: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
142: } else {
143: if (herm) {
144: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
145: MatSetOption(A,MAT_SPD,PETSC_TRUE);
146: } else {
147: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
148: MatSetOption(A,MAT_SPD,PETSC_FALSE);
149: }
150: MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
151: }
152: ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);
153: MatFactorSetSchurIS(F,is_schur);
154: ISDestroy(&is_schur);
155: if (use_lu) {
156: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
157: } else {
158: MatCholeskyFactorSymbolic(F,A,NULL,NULL);
159: }
161: for (nfact = 0; nfact < 3; nfact++) {
162: Mat AD;
164: if (!nfact) {
165: VecSetRandom(x,rand);
166: if (symm && herm) {
167: VecAbs(x);
168: }
169: MatDiagonalSet(A,x,ADD_VALUES);
170: }
171: if (use_lu) {
172: MatLUFactorNumeric(F,A,NULL);
173: } else {
174: MatCholeskyFactorNumeric(F,A,NULL);
175: }
176: MatFactorCreateSchurComplement(F,&S,NULL);
177: MatCreateVecs(S,&xschur,&bschur);
178: VecDuplicate(xschur,&uschur);
179: if (nfact == 1) {
180: MatFactorInvertSchurComplement(F);
181: }
182: for (nsolve = 0; nsolve < 2; nsolve++) {
183: VecSetRandom(x,rand);
184: VecCopy(x,u);
186: if (nsolve) {
187: MatMult(A,x,b);
188: MatSolve(F,b,x);
189: } else {
190: MatMultTranspose(A,x,b);
191: MatSolveTranspose(F,b,x);
192: }
193: /* Check the error */
194: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
195: VecNorm(u,NORM_2,&norm);
196: if (norm > tol) {
197: PetscReal resi;
198: if (nsolve) {
199: MatMult(A,x,u); /* u = A*x */
200: } else {
201: MatMultTranspose(A,x,u); /* u = A*x */
202: }
203: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
204: VecNorm(u,NORM_2,&resi);
205: if (nsolve) {
206: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
207: } else {
208: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
209: }
210: }
211: VecSetRandom(xschur,rand);
212: VecCopy(xschur,uschur);
213: if (nsolve) {
214: MatMult(S,xschur,bschur);
215: MatFactorSolveSchurComplement(F,bschur,xschur);
216: } else {
217: MatMultTranspose(S,xschur,bschur);
218: MatFactorSolveSchurComplementTranspose(F,bschur,xschur);
219: }
220: /* Check the error */
221: VecAXPY(uschur,-1.0,xschur); /* u <- (-1.0)x + u */
222: VecNorm(uschur,NORM_2,&norm);
223: if (norm > tol) {
224: PetscReal resi;
225: if (nsolve) {
226: MatMult(S,xschur,uschur); /* u = A*x */
227: } else {
228: MatMultTranspose(S,xschur,uschur); /* u = A*x */
229: }
230: VecAXPY(uschur,-1.0,bschur); /* u <- (-1.0)b + u */
231: VecNorm(uschur,NORM_2,&resi);
232: if (nsolve) {
233: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
234: } else {
235: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
236: }
237: }
238: }
239: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);
240: if (!nfact) {
241: MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);
242: } else {
243: MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);
244: }
245: MatDestroy(&AD);
246: for (nsolve = 0; nsolve < 2; nsolve++) {
247: MatMatSolve(F,RHS,X);
249: /* Check the error */
250: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
251: MatNorm(X,NORM_FROBENIUS,&norm);
252: if (norm > tol) {
253: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
254: }
255: }
256: if (isolver == 0) {
257: Mat spRHS,spRHST,RHST;
259: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
260: MatConvert(RHST,MATSEQAIJ,MAT_INITIAL_MATRIX,&spRHST);
261: MatCreateTranspose(spRHST,&spRHS);
262: for (nsolve = 0; nsolve < 2; nsolve++) {
263: MatMatSolve(F,spRHS,X);
265: /* Check the error */
266: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
267: MatNorm(X,NORM_FROBENIUS,&norm);
268: if (norm > tol) {
269: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) sparse MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
270: }
271: }
272: MatDestroy(&spRHST);
273: MatDestroy(&spRHS);
274: MatDestroy(&RHST);
275: }
276: MatDestroy(&S);
277: VecDestroy(&xschur);
278: VecDestroy(&bschur);
279: VecDestroy(&uschur);
280: }
281: /* Free data structures */
282: MatDestroy(&A);
283: MatDestroy(&C);
284: MatDestroy(&F);
285: MatDestroy(&X);
286: MatDestroy(&RHS);
287: PetscRandomDestroy(&rand);
288: VecDestroy(&x);
289: VecDestroy(&b);
290: VecDestroy(&u);
291: PetscFinalize();
292: return ierr;
293: }
296: /*TEST
298: test:
299: suffix: mkl_pardiso
300: requires: mkl_pardiso double !complex !define(PETSC_USE_64BIT_INDICES)
301: args: -solver 1
303: test:
304: suffix: mkl_pardiso_1
305: requires: mkl_pardiso double !complex !define(PETSC_USE_64BIT_INDICES)
306: args: -symmetric_solve -solver 1
308: test:
309: suffix: mkl_pardiso_3
310: requires: mkl_pardiso double !complex !define(PETSC_USE_64BIT_INDICES)
311: args: -symmetric_solve -hermitian_solve -solver 1
313: test:
314: suffix: mumps
315: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
316: args: -solver 0
318: test:
319: suffix: mumps_2
320: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
321: args: -symmetric_solve -solver 0
323: test:
324: suffix: mumps_3
325: requires: mumps double !complex !define(PETSC_USE_64BIT_INDICES)
326: args: -symmetric_solve -hermitian_solve -solver 0
328: TEST*/