Actual source code: ex192.c
petsc-3.8.4 2018-03-24
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: sprintf(file,PETSC_DIR);
38: PetscStrcat(file,"/share/petsc/datafiles/matrices/");
39: if (symm) {
40: #if defined (PETSC_USE_COMPLEX)
41: PetscStrcat(file,"hpd-complex-");
42: #else
43: PetscStrcat(file,"spd-real-");
44: #endif
45: } else {
46: #if defined (PETSC_USE_COMPLEX)
47: PetscStrcat(file,"nh-complex-");
48: #else
49: PetscStrcat(file,"ns-real-");
50: #endif
51: }
52: #if defined(PETSC_USE_64BIT_INDICES)
53: PetscStrcat(file,"int64-");
54: #else
55: PetscStrcat(file,"int32-");
56: #endif
57: #if defined (PETSC_USE_REAL_SINGLE)
58: PetscStrcat(file,"float32");
59: #else
60: PetscStrcat(file,"float64");
61: #endif
62: }
63: /* Load matrix A */
64: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
65: MatCreate(PETSC_COMM_WORLD,&A);
66: MatLoad(A,fd);
67: PetscViewerDestroy(&fd);
68: MatGetSize(A,&m,&n);
69: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
71: /* Create dense matrix C and X; C holds true solution with identical colums */
72: nrhs = 2;
73: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
74: MatCreate(PETSC_COMM_WORLD,&C);
75: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
76: MatSetType(C,MATDENSE);
77: MatSetFromOptions(C);
78: MatSetUp(C);
80: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
81: PetscRandomSetFromOptions(rand);
82: MatSetRandom(C,rand);
83: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
85: /* Create vectors */
86: VecCreate(PETSC_COMM_WORLD,&x);
87: VecSetSizes(x,n,PETSC_DECIDE);
88: VecSetFromOptions(x);
89: VecDuplicate(x,&b);
90: VecDuplicate(x,&u); /* save the true solution */
92: PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);
93: switch (isolver) {
94: #if defined(PETSC_HAVE_MUMPS)
95: case 0:
96: PetscStrcpy(solver,MATSOLVERMUMPS);
97: break;
98: #endif
99: #if defined(PETSC_HAVE_MKL_PARDISO)
100: case 1:
101: PetscStrcpy(solver,MATSOLVERMKL_PARDISO);
102: break;
103: #endif
104: default:
105: PetscStrcpy(solver,MATSOLVERPETSC);
106: break;
107: }
109: #if defined (PETSC_USE_COMPLEX)
110: if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
111: PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
112: PetscScalar val = -1.0;
113: val = val + im;
114: MatSetValue(A,1,0,val,INSERT_VALUES);
115: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
117: }
118: #endif
120: PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);
121: if (sratio < 0. || sratio > 1.) {
122: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
123: }
124: size_schur = (PetscInt)(sratio*m);
126: 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);
128: /* Test LU/Cholesky Factorization */
129: use_lu = PETSC_FALSE;
130: if (!symm) use_lu = PETSC_TRUE;
131: #if defined (PETSC_USE_COMPLEX)
132: if (isolver == 1) use_lu = PETSC_TRUE;
133: #endif
135: if (herm && !use_lu) { /* test also conversion routines inside the solver packages */
136: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
137: MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);
138: }
141: if (use_lu) {
142: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
143: } else {
144: if (herm) {
145: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
146: MatSetOption(A,MAT_SPD,PETSC_TRUE);
147: } else {
148: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
149: MatSetOption(A,MAT_SPD,PETSC_FALSE);
150: }
151: MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
152: }
153: ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);
154: MatFactorSetSchurIS(F,is_schur);
155: ISDestroy(&is_schur);
156: if (use_lu) {
157: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
158: } else {
159: MatCholeskyFactorSymbolic(F,A,NULL,NULL);
160: }
162: for (nfact = 0; nfact < 3; nfact++) {
163: Mat AD;
165: if (!nfact) {
166: VecSetRandom(x,rand);
167: if (symm && herm) {
168: VecAbs(x);
169: }
170: MatDiagonalSet(A,x,ADD_VALUES);
171: }
172: if (use_lu) {
173: MatLUFactorNumeric(F,A,NULL);
174: } else {
175: MatCholeskyFactorNumeric(F,A,NULL);
176: }
177: MatFactorCreateSchurComplement(F,&S,NULL);
178: MatCreateVecs(S,&xschur,&bschur);
179: VecDuplicate(xschur,&uschur);
180: if (nfact == 1) {
181: MatFactorInvertSchurComplement(F);
182: }
183: for (nsolve = 0; nsolve < 2; nsolve++) {
184: VecSetRandom(x,rand);
185: VecCopy(x,u);
187: if (nsolve) {
188: MatMult(A,x,b);
189: MatSolve(F,b,x);
190: } else {
191: MatMultTranspose(A,x,b);
192: MatSolveTranspose(F,b,x);
193: }
194: /* Check the error */
195: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
196: VecNorm(u,NORM_2,&norm);
197: if (norm > tol) {
198: PetscReal resi;
199: if (nsolve) {
200: MatMult(A,x,u); /* u = A*x */
201: } else {
202: MatMultTranspose(A,x,u); /* u = A*x */
203: }
204: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
205: VecNorm(u,NORM_2,&resi);
206: if (nsolve) {
207: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
208: } else {
209: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
210: }
211: }
212: VecSetRandom(xschur,rand);
213: VecCopy(xschur,uschur);
214: if (nsolve) {
215: MatMult(S,xschur,bschur);
216: MatFactorSolveSchurComplement(F,bschur,xschur);
217: } else {
218: MatMultTranspose(S,xschur,bschur);
219: MatFactorSolveSchurComplementTranspose(F,bschur,xschur);
220: }
221: /* Check the error */
222: VecAXPY(uschur,-1.0,xschur); /* u <- (-1.0)x + u */
223: VecNorm(uschur,NORM_2,&norm);
224: if (norm > tol) {
225: PetscReal resi;
226: if (nsolve) {
227: MatMult(S,xschur,uschur); /* u = A*x */
228: } else {
229: MatMultTranspose(S,xschur,uschur); /* u = A*x */
230: }
231: VecAXPY(uschur,-1.0,bschur); /* u <- (-1.0)b + u */
232: VecNorm(uschur,NORM_2,&resi);
233: if (nsolve) {
234: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
235: } else {
236: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
237: }
238: }
239: }
240: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);
241: if (!nfact) {
242: MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);
243: } else {
244: MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);
245: }
246: MatDestroy(&AD);
247: for (nsolve = 0; nsolve < 2; nsolve++) {
248: MatMatSolve(F,RHS,X);
250: /* Check the error */
251: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
252: MatNorm(X,NORM_FROBENIUS,&norm);
253: if (norm > tol) {
254: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
255: }
256: }
257: if (isolver == 0) {
258: Mat spRHS,spRHST,RHST;
260: MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);
261: MatConvert(RHST,MATSEQAIJ,MAT_INITIAL_MATRIX,&spRHST);
262: MatCreateTranspose(spRHST,&spRHS);
263: for (nsolve = 0; nsolve < 2; nsolve++) {
264: MatMatSolve(F,spRHS,X);
266: /* Check the error */
267: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
268: MatNorm(X,NORM_FROBENIUS,&norm);
269: if (norm > tol) {
270: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) sparse MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
271: }
272: }
273: MatDestroy(&spRHST);
274: MatDestroy(&spRHS);
275: MatDestroy(&RHST);
276: }
277: MatDestroy(&S);
278: VecDestroy(&xschur);
279: VecDestroy(&bschur);
280: VecDestroy(&uschur);
281: }
282: /* Free data structures */
283: MatDestroy(&A);
284: MatDestroy(&C);
285: MatDestroy(&F);
286: MatDestroy(&X);
287: MatDestroy(&RHS);
288: PetscRandomDestroy(&rand);
289: VecDestroy(&x);
290: VecDestroy(&b);
291: VecDestroy(&u);
292: PetscFinalize();
293: return ierr;
294: }