Actual source code: ex192.c
petsc-3.6.4 2016-04-12
2: static char help[] = "Tests MatSolve() and MatMatSolve() with mumps sequential solver in Schur complement mode.\n\
3: Example: mpiexec -n 1 ./ex192 -f <matrix binary file> -nrhs 4 -symmetric_solve -hermitian_solve -sratio 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;
15: PetscMPIInt size;
16: PetscInt icntl19,size_schur,*idxs_schur,i,m,n,nfact,nsolve,nrhs;
17: PetscReal norm,tol=1.e-12;
18: PetscRandom rand;
19: PetscBool flg,herm,symm;
20: PetscReal sratio = 5.1/12.;
21: PetscViewer fd; /* viewer */
22: char file[PETSC_MAX_PATH_LEN]; /* input file name */
24: PetscInitialize(&argc,&args,(char*)0,help);
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,"-symmetric_solve",&symm,NULL);
31: PetscOptionsGetBool(NULL,"-hermitian_solve",&herm,NULL);
32: if (herm) symm = PETSC_TRUE;
34: /* Determine file from which we read the matrix A */
35: PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
36: if (!flg) { /* 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);
68: #if defined (PETSC_USE_COMPLEX)
69: if (symm & !flg) { /* MUMPS (5.0.0) does not have support for hermitian matrices, so make them symmetric */
70: PetscScalar im = PetscSqrtScalar((PetscScalar)-1.);
71: PetscScalar val = -1.0;
72: val = val + im;
73: MatSetValue(A,1,0,val,INSERT_VALUES);
74: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
76: }
77: #endif
78: PetscViewerDestroy(&fd);
79: MatGetSize(A,&m,&n);
80: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
82: /* Create dense matrix C and X; C holds true solution with identical colums */
83: nrhs = 2;
84: PetscOptionsGetInt(NULL,"-nrhs",&nrhs,NULL);
85: MatCreate(PETSC_COMM_WORLD,&C);
86: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
87: MatSetType(C,MATDENSE);
88: MatSetFromOptions(C);
89: MatSetUp(C);
91: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
92: PetscRandomSetFromOptions(rand);
93: MatSetRandom(C,rand);
94: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
96: /* Create vectors */
97: VecCreate(PETSC_COMM_WORLD,&x);
98: VecSetSizes(x,n,PETSC_DECIDE);
99: VecSetFromOptions(x);
100: VecDuplicate(x,&b);
101: VecDuplicate(x,&u); /* save the true solution */
103: /* Test LU/Cholesky Factorization */
104: if (!symm) {
105: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
106: } else {
107: if (herm) {
108: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
109: MatSetOption(A,MAT_SPD,PETSC_TRUE);
110: } else {
111: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
112: MatSetOption(A,MAT_SPD,PETSC_FALSE);
113: }
114: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);
115: }
116: PetscOptionsGetReal(NULL,"-schur_ratio",&sratio,NULL);
117: if (sratio < 0. || sratio > 1.) {
118: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid ratio for schur degrees of freedom %f", sratio);
119: }
120: size_schur = (PetscInt)(sratio*m);
121: PetscMalloc1(size_schur,&idxs_schur);
122: for (i=0;i<size_schur;i++) {
123: idxs_schur[i] = m-size_schur+i+1; /* fortran like */
124: }
125: MatMumpsSetSchurIndices(F,size_schur,idxs_schur);
126: PetscFree(idxs_schur);
127: if (!symm) {
128: MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
129: } else {
130: MatCholeskyFactorSymbolic(F,A,NULL,NULL);
131: }
132: MatMumpsGetIcntl(F,19,&icntl19);
133: PetscPrintf(PETSC_COMM_SELF,"Solving: nrhs %d, sym %d, herm %d, size schur %d, size mat %d, icntl19 %d\n",nrhs,symm,herm,size_schur,m,icntl19);
135: for (nfact = 0; nfact < 3; nfact++) {
137: if (!nfact) {
138: VecSetRandom(x,rand);
139: if (symm && herm) {
140: VecAbs(x);
141: }
142: MatDiagonalSet(A,x,ADD_VALUES);
143: }
144: if (!symm) {
145: MatLUFactorNumeric(F,A,NULL);
146: } else {
147: MatCholeskyFactorNumeric(F,A,NULL);
148: }
149: if (icntl19) {
150: MatMumpsCreateSchurComplement(F,&S);
151: MatCreateVecs(S,&xschur,&bschur);
152: VecDuplicate(xschur,&uschur);
153: if (nfact == 1) {
154: MatMumpsInvertSchurComplement(F);
155: }
156: }
157: for (nsolve = 0; nsolve < 2; nsolve++) {
158: VecSetRandom(x,rand);
159: VecCopy(x,u);
161: if (nsolve) {
162: MatMult(A,x,b);
163: MatSolve(F,b,x);
164: } else {
165: MatMultTranspose(A,x,b);
166: MatSolveTranspose(F,b,x);
167: }
168: /* Check the error */
169: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
170: VecNorm(u,NORM_2,&norm);
171: if (norm > tol) {
172: PetscReal resi;
173: if (nsolve) {
174: MatMult(A,x,u); /* u = A*x */
175: } else {
176: MatMultTranspose(A,x,u); /* u = A*x */
177: }
178: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
179: VecNorm(u,NORM_2,&resi);
180: if (nsolve) {
181: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolve: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
182: } else {
183: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatSolveTranspose: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
184: }
185: }
186: if (icntl19) {
187: VecSetRandom(xschur,rand);
188: VecCopy(xschur,uschur);
189: if (nsolve) {
190: MatMult(S,xschur,bschur);
191: MatMumpsSolveSchurComplement(F,bschur,xschur);
192: } else {
193: MatMultTranspose(S,xschur,bschur);
194: MatMumpsSolveSchurComplementTranspose(F,bschur,xschur);
195: }
196: /* Check the error */
197: VecAXPY(uschur,-1.0,xschur); /* u <- (-1.0)x + u */
198: VecNorm(uschur,NORM_2,&norm);
199: if (norm > tol) {
200: PetscReal resi;
201: if (nsolve) {
202: MatMult(S,xschur,uschur); /* u = A*x */
203: } else {
204: MatMultTranspose(S,xschur,uschur); /* u = A*x */
205: }
206: VecAXPY(uschur,-1.0,bschur); /* u <- (-1.0)b + u */
207: VecNorm(uschur,NORM_2,&resi);
208: if (nsolve) {
209: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatMumpsSolveSchurComplement: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
210: } else {
211: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatMumpsSolveSchurComplementTranspose: Norm of error %g, residual %f\n",nfact,nsolve,norm,resi);
212: }
213: }
214: }
215: }
217: if (!nfact) {
218: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
219: } else {
220: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
221: }
222: for (nsolve = 0; nsolve < 2; nsolve++) {
223: MatMatSolve(F,RHS,X);
225: /* Check the error */
226: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
227: MatNorm(X,NORM_FROBENIUS,&norm);
228: if (norm > tol) {
229: PetscPrintf(PETSC_COMM_SELF,"(f %d, s %d) MatMatSolve: Norm of error %g\n",nfact,nsolve,norm);
230: }
231: }
232: if (icntl19) {
233: MatDestroy(&S);
234: VecDestroy(&xschur);
235: VecDestroy(&bschur);
236: VecDestroy(&uschur);
237: }
238: }
239: /* Free data structures */
240: MatDestroy(&A);
241: MatDestroy(&C);
242: MatDestroy(&F);
243: MatDestroy(&X);
244: MatDestroy(&RHS);
245: PetscRandomDestroy(&rand);
246: VecDestroy(&x);
247: VecDestroy(&b);
248: VecDestroy(&u);
249: PetscFinalize();
250: return 0;
251: }