Actual source code: ex192.c
petsc-3.7.3 2016-08-01
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>
9: int main(int argc,char **args)
10: {
11: Mat A,RHS,C,F,X,S;
12: Vec u,x,b;
13: Vec xschur,bschur,uschur;
14: IS is_schur;
16: PetscMPIInt size;
17: PetscInt isolver=0,size_schur,m,n,nfact,nsolve,nrhs;
18: PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON;
19: PetscRandom rand;
20: PetscBool data_provided,herm,symm,use_lu;
21: PetscReal sratio = 5.1/12.;
22: PetscViewer fd; /* viewer */
23: char solver[256];
24: char file[PETSC_MAX_PATH_LEN]; /* input file name */
26: PetscInitialize(&argc,&args,(char*)0,help);
27: MPI_Comm_size(PETSC_COMM_WORLD, &size);
28: if (size > 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor test");
29: /* Determine which type of solver we want to test for */
30: herm = PETSC_FALSE;
31: symm = PETSC_FALSE;
32: PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);
33: PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);
34: if (herm) symm = PETSC_TRUE;
36: /* Determine file from which we read the matrix A */
37: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&data_provided);
38: if (!data_provided) { /* get matrices from PETSc distribution */
39: sprintf(file,PETSC_DIR);
40: PetscStrcat(file,"/share/petsc/datafiles/matrices/");
41: if (symm) {
42: #if defined (PETSC_USE_COMPLEX)
43: PetscStrcat(file,"hpd-complex-");
44: #else
45: PetscStrcat(file,"spd-real-");
46: #endif
47: } else {
48: #if defined (PETSC_USE_COMPLEX)
49: PetscStrcat(file,"nh-complex-");
50: #else
51: PetscStrcat(file,"ns-real-");
52: #endif
53: }
54: #if defined(PETSC_USE_64BIT_INDICES)
55: PetscStrcat(file,"int64-");
56: #else
57: PetscStrcat(file,"int32-");
58: #endif
59: #if defined (PETSC_USE_REAL_SINGLE)
60: PetscStrcat(file,"float32");
61: #else
62: PetscStrcat(file,"float64");
63: #endif
64: }
65: /* Load matrix A */
66: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
67: MatCreate(PETSC_COMM_WORLD,&A);
68: MatLoad(A,fd);
69: PetscViewerDestroy(&fd);
70: MatGetSize(A,&m,&n);
71: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
73: /* Create dense matrix C and X; C holds true solution with identical colums */
74: nrhs = 2;
75: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
76: MatCreate(PETSC_COMM_WORLD,&C);
77: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
78: MatSetType(C,MATDENSE);
79: MatSetFromOptions(C);
80: MatSetUp(C);
82: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
83: PetscRandomSetFromOptions(rand);
84: MatSetRandom(C,rand);
85: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
87: /* Create vectors */
88: VecCreate(PETSC_COMM_WORLD,&x);
89: VecSetSizes(x,n,PETSC_DECIDE);
90: VecSetFromOptions(x);
91: VecDuplicate(x,&b);
92: VecDuplicate(x,&u); /* save the true solution */
94: PetscOptionsGetInt(NULL,NULL,"-solver",&isolver,NULL);
95: switch (isolver) {
96: #if defined(PETSC_HAVE_MUMPS)
97: case 0:
98: PetscStrcpy(solver,MATSOLVERMUMPS);
99: break;
100: #endif
101: #if defined(PETSC_HAVE_MKL_PARDISO)
102: case 1:
103: PetscStrcpy(solver,MATSOLVERMKL_PARDISO);
104: break;
105: #endif
106: default:
107: PetscStrcpy(solver,MATSOLVERPETSC);
108: break;
109: }
111: #if defined (PETSC_USE_COMPLEX)
112: if (isolver == 0 && symm && !data_provided) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
113: PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
114: PetscScalar val = -1.0;
115: val = val + im;
116: MatSetValue(A,1,0,val,INSERT_VALUES);
117: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
118: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
119: }
120: #endif
122: PetscOptionsGetReal(NULL,NULL,"-schur_ratio",&sratio,NULL);
123: if (sratio < 0. || sratio > 1.) {
124: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
125: }
126: size_schur = (PetscInt)(sratio*m);
128: 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);
130: /* Test LU/Cholesky Factorization */
131: use_lu = PETSC_FALSE;
132: if (!symm) use_lu = PETSC_TRUE;
133: #if defined (PETSC_USE_COMPLEX)
134: if (isolver == 1) use_lu = PETSC_TRUE;
135: #endif
137: if (herm && !use_lu) { /* test also conversion routines inside the solver packages */
138: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
139: MatConvert(A,MATSEQSBAIJ,MAT_INPLACE_MATRIX,&A);
140: }
143: if (use_lu) {
144: MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
145: } else {
146: if (herm) {
147: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
148: MatSetOption(A,MAT_SPD,PETSC_TRUE);
149: } else {
150: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
151: MatSetOption(A,MAT_SPD,PETSC_FALSE);
152: }
153: MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
154: }
155: ISCreateStride(PETSC_COMM_SELF,size_schur,m-size_schur,1,&is_schur);
156: MatFactorSetSchurIS(F,is_schur);
157: ISDestroy(&is_schur);
158: if (use_lu) {
159: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
160: } else {
161: MatCholeskyFactorSymbolic(F,A,NULL,NULL);
162: }
164: for (nfact = 0; nfact < 3; nfact++) {
165: Mat AD;
167: if (!nfact) {
168: VecSetRandom(x,rand);
169: if (symm && herm) {
170: VecAbs(x);
171: }
172: MatDiagonalSet(A,x,ADD_VALUES);
173: }
174: if (use_lu) {
175: MatLUFactorNumeric(F,A,NULL);
176: } else {
177: MatCholeskyFactorNumeric(F,A,NULL);
178: }
179: MatFactorCreateSchurComplement(F,&S);
180: MatCreateVecs(S,&xschur,&bschur);
181: VecDuplicate(xschur,&uschur);
182: if (nfact == 1) {
183: MatFactorInvertSchurComplement(F);
184: }
185: for (nsolve = 0; nsolve < 2; nsolve++) {
186: VecSetRandom(x,rand);
187: VecCopy(x,u);
189: if (nsolve) {
190: MatMult(A,x,b);
191: MatSolve(F,b,x);
192: } else {
193: MatMultTranspose(A,x,b);
194: MatSolveTranspose(F,b,x);
195: }
196: /* Check the error */
197: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
198: VecNorm(u,NORM_2,&norm);
199: if (norm > tol) {
200: PetscReal resi;
201: if (nsolve) {
202: MatMult(A,x,u); /* u = A*x */
203: } else {
204: MatMultTranspose(A,x,u); /* u = A*x */
205: }
206: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
207: VecNorm(u,NORM_2,&resi);
208: if (nsolve) {
209: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolve error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
210: } else {
211: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolveTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
212: }
213: }
214: VecSetRandom(xschur,rand);
215: VecCopy(xschur,uschur);
216: if (nsolve) {
217: MatMult(S,xschur,bschur);
218: MatFactorSolveSchurComplement(F,bschur,xschur);
219: } else {
220: MatMultTranspose(S,xschur,bschur);
221: MatFactorSolveSchurComplementTranspose(F,bschur,xschur);
222: }
223: /* Check the error */
224: VecAXPY(uschur,-1.0,xschur); /* u <- (-1.0)x + u */
225: VecNorm(uschur,NORM_2,&norm);
226: if (norm > tol) {
227: PetscReal resi;
228: if (nsolve) {
229: MatMult(S,xschur,uschur); /* u = A*x */
230: } else {
231: MatMultTranspose(S,xschur,uschur); /* u = A*x */
232: }
233: VecAXPY(uschur,-1.0,bschur); /* u <- (-1.0)b + u */
234: VecNorm(uschur,NORM_2,&resi);
235: if (nsolve) {
236: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatFactorSolveSchurComplement error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
237: } else {
238: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatFactorSolveSchurComplementTranspose error: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
239: }
240: }
241: }
242: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&AD);
243: if (!nfact) {
244: MatMatMult(AD,C,MAT_INITIAL_MATRIX,2.0,&RHS);
245: } else {
246: MatMatMult(AD,C,MAT_REUSE_MATRIX,2.0,&RHS);
247: }
248: MatDestroy(&AD);
249: for (nsolve = 0; nsolve < 2; nsolve++) {
250: MatMatSolve(F,RHS,X);
252: /* Check the error */
253: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
254: MatNorm(X,NORM_FROBENIUS,&norm);
255: if (norm > tol) {
256: PetscPrintf(PETSC_COMM_SELF,"(f %D, s %D) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
257: }
258: }
259: MatDestroy(&S);
260: VecDestroy(&xschur);
261: VecDestroy(&bschur);
262: VecDestroy(&uschur);
263: }
264: /* Free data structures */
265: MatDestroy(&A);
266: MatDestroy(&C);
267: MatDestroy(&F);
268: MatDestroy(&X);
269: MatDestroy(&RHS);
270: PetscRandomDestroy(&rand);
271: VecDestroy(&x);
272: VecDestroy(&b);
273: VecDestroy(&u);
274: PetscFinalize();
275: return 0;
276: }