Actual source code: ex125.c
petsc-3.7.3 2016-08-01
2: static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
3: Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n";
5: #include <petscmat.h>
9: int main(int argc,char **args)
10: {
11: Mat A,RHS,C,F,X;
12: Vec u,x,b;
14: PetscMPIInt rank,size;
15: PetscInt i,m,n,nfact,nsolve,nrhs,ipack=0;
16: PetscScalar *array,rval;
17: PetscReal norm,tol=1.e-12;
18: IS perm,iperm;
19: MatFactorInfo info;
20: PetscRandom rand;
21: PetscBool flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE;
22: PetscViewer fd; /* viewer */
23: char file[PETSC_MAX_PATH_LEN]; /* input file name */
25: PetscInitialize(&argc,&args,(char*)0,help);
26: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
27: MPI_Comm_size(PETSC_COMM_WORLD, &size);
29: /* Determine file from which we read the matrix A */
30: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
31: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
33: /* Load matrix A */
34: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
35: MatCreate(PETSC_COMM_WORLD,&A);
36: MatLoad(A,fd);
37: PetscViewerDestroy(&fd);
38: MatGetLocalSize(A,&m,&n);
39: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
41: /* Create dense matrix C and X; C holds true solution with identical colums */
42: nrhs = 2;
43: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
44: if (!rank) printf("ex125: nrhs %d\n",nrhs);
45: MatCreate(PETSC_COMM_WORLD,&C);
46: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
47: MatSetType(C,MATDENSE);
48: MatSetFromOptions(C);
49: MatSetUp(C);
50:
51: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
52: PetscRandomSetFromOptions(rand);
53: /* #define DEBUGEX */
54: #if defined(DEBUGEX)
55: {
56: PetscInt row,j,M,cols[nrhs];
57: PetscScalar vals[nrhs];
58: MatGetSize(A,&M,NULL);
59: if (!rank) {
60: for (j=0; j<nrhs; j++) cols[j] = j;
61: for (row = 0; row < M; row++){
62: for (j=0; j<nrhs; j++) vals[j] = row;
63: MatSetValues(C,1,&row,nrhs,cols,vals,INSERT_VALUES);
64: }
65: }
66: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
68: }
69: #else
70: MatSetRandom(C,rand);
71: #endif
72: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
74: /* Create vectors */
75: VecCreate(PETSC_COMM_WORLD,&x);
76: VecSetSizes(x,n,PETSC_DECIDE);
77: VecSetFromOptions(x);
78: VecDuplicate(x,&b);
79: VecDuplicate(x,&u); /* save the true solution */
81: /* Test LU Factorization */
82: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
83: /*ISView(perm,PETSC_VIEWER_STDOUT_WORLD);*/
84: /*ISView(perm,PETSC_VIEWER_STDOUT_SELF);*/
86: PetscOptionsGetInt(NULL,NULL,"-mat_solver_package",&ipack,NULL);
87: switch (ipack) {
88: #if defined(PETSC_HAVE_SUPERLU)
89: case 0:
90: if (!rank) printf(" SUPERLU LU:\n");
91: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
92: break;
93: #endif
94: #if defined(PETSC_HAVE_SUPERLU_DIST)
95: case 1:
96: if (!rank) printf(" SUPERLU_DIST LU:\n");
97: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
98: break;
99: #endif
100: #if defined(PETSC_HAVE_MUMPS)
101: case 2:
102: if (!rank) printf(" MUMPS LU:\n");
103: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
104: {
105: /* test mumps options */
106: PetscInt icntl;
107: PetscReal cntl;
108:
109: icntl = 2; /* sequential matrix ordering */
110: MatMumpsSetIcntl(F,7,icntl);
112: cntl = 1.e-6; /* threshhold for row pivot detection */
113: MatMumpsSetIcntl(F,24,1);
114: MatMumpsSetCntl(F,3,cntl);
115: }
116: break;
117: #endif
118: #if defined(PETSC_HAVE_MKL_PARDISO)
119: case 3:
120: if (!rank) printf(" MKL_PARDISO LU:\n");
121: MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);
122: break;
123: #endif
124: default:
125: if (!rank) printf(" PETSC LU:\n");
126: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
127: }
129: MatFactorInfoInitialize(&info);
130: info.fill = 5.0;
131: info.shifttype = (PetscReal) MAT_SHIFT_NONE;
132: MatLUFactorSymbolic(F,A,perm,iperm,&info);
134: for (nfact = 0; nfact < 2; nfact++) {
135: if (!rank) printf(" %d-the LU numfactorization \n",nfact);
136: MatLUFactorNumeric(F,A,&info);
138: /* Test MatMatSolve() */
139: if (testMatMatSolve) {
140: if (!nfact) {
141: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
142: } else {
143: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
144: }
145: for (nsolve = 0; nsolve < 2; nsolve++) {
146: if (!rank) printf(" %d-the MatMatSolve \n",nsolve);
147: MatMatSolve(F,RHS,X);
149: /* Check the error */
150: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
151: MatNorm(X,NORM_FROBENIUS,&norm);
152: if (norm > tol) {
153: if (!rank) {
154: PetscPrintf(PETSC_COMM_SELF,"%d-the MatMatSolve: Norm of error %g, nsolve %d\n",nsolve,norm,nsolve);
155: }
156: }
157: }
158: }
160: /* Test MatSolve() */
161: if (testMatSolve) {
162: for (nsolve = 0; nsolve < 2; nsolve++) {
163: VecGetArray(x,&array);
164: for (i=0; i<m; i++) {
165: PetscRandomGetValue(rand,&rval);
166: array[i] = rval;
167: }
168: VecRestoreArray(x,&array);
169: VecCopy(x,u);
170: MatMult(A,x,b);
172: if (!rank) printf(" %d-the MatSolve \n",nsolve);
173: MatSolve(F,b,x);
175: /* Check the error */
176: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
177: VecNorm(u,NORM_2,&norm);
178: if (norm > tol) {
179: PetscReal resi;
180: MatMult(A,x,u); /* u = A*x */
181: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
182: VecNorm(u,NORM_2,&resi);
183: if (!rank) {
184: PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g, resi %g, LU numfact %d\n",norm,resi,nfact);
185: }
186: }
187: }
188: }
189: }
191: /* Free data structures */
192: MatDestroy(&A);
193: MatDestroy(&C);
194: MatDestroy(&F);
195: MatDestroy(&X);
196: if (testMatMatSolve) {
197: MatDestroy(&RHS);
198: }
200: PetscRandomDestroy(&rand);
201: ISDestroy(&perm);
202: ISDestroy(&iperm);
203: VecDestroy(&x);
204: VecDestroy(&b);
205: VecDestroy(&u);
206: PetscFinalize();
207: return 0;
208: }