xref: /petsc/src/mat/tests/ex125.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce) 
1 static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
2 Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n";
3 
4 #include <petscmat.h>
5 
6 int main(int argc,char **args)
7 {
8   Mat            A,RHS,C,F,X;
9   Vec            u,x,b;
10   PetscErrorCode ierr;
11   PetscMPIInt    size;
12   PetscInt       m,n,nfact,nsolve,nrhs,ipack=0;
13   PetscReal      norm,tol=1.e-10;
14   IS             perm,iperm;
15   MatFactorInfo  info;
16   PetscRandom    rand;
17   PetscBool      flg,testMatSolve=PETSC_TRUE,testMatMatSolve=PETSC_TRUE;
18   PetscBool      chol=PETSC_FALSE,view=PETSC_FALSE,matsolvexx = PETSC_FALSE;
19 #if defined(PETSC_HAVE_MUMPS)
20   PetscBool      test_mumps_opts=PETSC_FALSE;
21 #endif
22   PetscViewer    fd;              /* viewer */
23   char           file[PETSC_MAX_PATH_LEN]; /* input file name */
24 
25   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
26   ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr);
27 
28   /* Determine file from which we read the matrix A */
29   ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);CHKERRQ(ierr);
30   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
31 
32   /* Load matrix A */
33   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
34   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
35   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
36   ierr = MatLoad(A,fd);CHKERRQ(ierr);
37   ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
38   ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
39   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%D, %D)", m, n);
40 
41   /* if A is symmetric, set its flag -- required by MatGetInertia() */
42   ierr = MatIsSymmetric(A,0.0,&flg);CHKERRQ(ierr);
43 
44   ierr = MatViewFromOptions(A,NULL,"-A_view");CHKERRQ(ierr);
45 
46   /* Create dense matrix C and X; C holds true solution with identical colums */
47   nrhs = 2;
48   ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr);
49   ierr = PetscPrintf(PETSC_COMM_WORLD,"ex125: nrhs %D\n",nrhs);CHKERRQ(ierr);
50   ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
51   ierr = MatSetOptionsPrefix(C,"rhs_");CHKERRQ(ierr);
52   ierr = MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);CHKERRQ(ierr);
53   ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr);
54   ierr = MatSetFromOptions(C);CHKERRQ(ierr);
55   ierr = MatSetUp(C);CHKERRQ(ierr);
56 
57   ierr = PetscOptionsGetBool(NULL,NULL,"-view_factor",&view,NULL);CHKERRQ(ierr);
58   ierr = PetscOptionsGetBool(NULL,NULL,"-test_matmatsolve",&testMatMatSolve,NULL);CHKERRQ(ierr);
59   ierr = PetscOptionsGetBool(NULL,NULL,"-cholesky",&chol,NULL);CHKERRQ(ierr);
60 #if defined(PETSC_HAVE_MUMPS)
61   ierr = PetscOptionsGetBool(NULL,NULL,"-test_mumps_opts",&test_mumps_opts,NULL);CHKERRQ(ierr);
62 #endif
63 
64   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr);
65   ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr);
66   ierr = MatSetRandom(C,rand);CHKERRQ(ierr);
67   ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);CHKERRQ(ierr);
68 
69   /* Create vectors */
70   ierr = MatCreateVecs(A,&x,&b);CHKERRQ(ierr);
71   ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* save the true solution */
72 
73   /* Test Factorization */
74   ierr = MatGetOrdering(A,MATORDERINGND,&perm,&iperm);CHKERRQ(ierr);
75 
76   ierr = PetscOptionsGetInt(NULL,NULL,"-mat_solver_type",&ipack,NULL);CHKERRQ(ierr);
77   switch (ipack) {
78 #if defined(PETSC_HAVE_SUPERLU)
79   case 0:
80     if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!");
81     ierr = PetscPrintf(PETSC_COMM_WORLD," SUPERLU LU:\n");CHKERRQ(ierr);
82     ierr = MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
83     matsolvexx = PETSC_TRUE;
84     break;
85 #endif
86 #if defined(PETSC_HAVE_SUPERLU_DIST)
87   case 1:
88     if (chol) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"SuperLU does not provide Cholesky!");
89     ierr = PetscPrintf(PETSC_COMM_WORLD," SUPERLU_DIST LU:\n");CHKERRQ(ierr);
90     ierr = MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
91     matsolvexx = PETSC_TRUE;
92     break;
93 #endif
94 #if defined(PETSC_HAVE_MUMPS)
95   case 2:
96     if (chol) {
97       ierr = PetscPrintf(PETSC_COMM_WORLD," MUMPS CHOLESKY:\n");CHKERRQ(ierr);
98       ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
99     } else {
100       ierr = PetscPrintf(PETSC_COMM_WORLD," MUMPS LU:\n");CHKERRQ(ierr);
101       ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
102     }
103     matsolvexx = PETSC_TRUE;
104     if (test_mumps_opts) {
105       /* test mumps options */
106       PetscInt  icntl;
107       PetscReal cntl;
108 
109       icntl = 2;        /* sequential matrix ordering */
110       ierr  = MatMumpsSetIcntl(F,7,icntl);CHKERRQ(ierr);
111 
112       cntl = 1.e-6; /* threshold for row pivot detection */
113       ierr = MatMumpsSetIcntl(F,24,1);CHKERRQ(ierr);
114       ierr = MatMumpsSetCntl(F,3,cntl);CHKERRQ(ierr);
115     }
116     break;
117 #endif
118 #if defined(PETSC_HAVE_MKL_PARDISO)
119   case 3:
120     if (chol) {
121       ierr = PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO CHOLESKY:\n");CHKERRQ(ierr);
122       ierr = MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
123     } else {
124       ierr = PetscPrintf(PETSC_COMM_WORLD," MKL_PARDISO LU:\n");CHKERRQ(ierr);
125       ierr = MatGetFactor(A,MATSOLVERMKL_PARDISO,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
126     }
127     break;
128 #endif
129 #if defined(PETSC_HAVE_CUDA)
130   case 4:
131     if (chol) {
132       ierr = PetscPrintf(PETSC_COMM_WORLD," CUSPARSE CHOLESKY:\n");CHKERRQ(ierr);
133       ierr = MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
134     } else {
135       ierr = PetscPrintf(PETSC_COMM_WORLD," CUSPARSE LU:\n");CHKERRQ(ierr);
136       ierr = MatGetFactor(A,MATSOLVERCUSPARSE,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
137     }
138     break;
139 #endif
140   default:
141     if (chol) {
142       ierr = PetscPrintf(PETSC_COMM_WORLD," PETSC CHOLESKY:\n");CHKERRQ(ierr);
143       ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
144     } else {
145       ierr = PetscPrintf(PETSC_COMM_WORLD," PETSC LU:\n");CHKERRQ(ierr);
146       ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
147     }
148     matsolvexx = PETSC_TRUE;
149   }
150 
151   ierr           = MatFactorInfoInitialize(&info);CHKERRQ(ierr);
152   info.fill      = 5.0;
153   info.shifttype = (PetscReal) MAT_SHIFT_NONE;
154   if (chol) {
155     ierr = MatCholeskyFactorSymbolic(F,A,perm,&info);CHKERRQ(ierr);
156   } else {
157     ierr = MatLUFactorSymbolic(F,A,perm,iperm,&info);CHKERRQ(ierr);
158   }
159 
160   for (nfact = 0; nfact < 2; nfact++) {
161     if (chol) {
162       ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the CHOLESKY numfactorization \n",nfact);CHKERRQ(ierr);
163       ierr = MatCholeskyFactorNumeric(F,A,&info);CHKERRQ(ierr);
164     } else {
165       ierr = PetscPrintf(PETSC_COMM_WORLD," %D-the LU numfactorization \n",nfact);CHKERRQ(ierr);
166       ierr = MatLUFactorNumeric(F,A,&info);CHKERRQ(ierr);
167     }
168     if (view) {
169       ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
170       ierr = MatView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
171       ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
172       view = PETSC_FALSE;
173     }
174 
175 #if defined(PETSC_HAVE_SUPERLU_DIST)
176     if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
177        -- input: matrix factor F; output: main diagonal of matrix U on all processes */
178       PetscInt    M;
179       PetscScalar *diag;
180 #if !defined(PETSC_USE_COMPLEX)
181       PetscInt nneg,nzero,npos;
182 #endif
183 
184       ierr = MatGetSize(F,&M,NULL);CHKERRQ(ierr);
185       ierr = PetscMalloc1(M,&diag);CHKERRQ(ierr);
186       ierr = MatSuperluDistGetDiagU(F,diag);CHKERRQ(ierr);
187       ierr = PetscFree(diag);CHKERRQ(ierr);
188 
189 #if !defined(PETSC_USE_COMPLEX)
190       /* Test MatGetInertia() */
191       ierr = MatGetInertia(F,&nneg,&nzero,&npos);CHKERRQ(ierr);
192       ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);CHKERRQ(ierr);
193 #endif
194     }
195 #endif
196 
197     /* Test MatMatSolve() */
198     if (testMatMatSolve) {
199       if (!nfact) {
200         ierr = MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);CHKERRQ(ierr);
201       } else {
202         ierr = MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr);
203       }
204       for (nsolve = 0; nsolve < 2; nsolve++) {
205         ierr = PetscPrintf(PETSC_COMM_WORLD,"   %D-the MatMatSolve \n",nsolve);CHKERRQ(ierr);
206         ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr);
207 
208         /* Check the error */
209         ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
210         ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
211         if (norm > tol) {
212           ierr = PetscPrintf(PETSC_COMM_WORLD,"%D-the MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);CHKERRQ(ierr);
213         }
214       }
215       if (matsolvexx) {
216         /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
217         ierr = MatCopy(RHS,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
218         ierr = MatMatSolve(F,X,X);CHKERRQ(ierr);
219         /* Check the error */
220         ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
221         ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
222         if (norm > tol) {
223           ierr = PetscPrintf(PETSC_COMM_WORLD,"MatMatSolve(F,RHS,RHS): Norm of error %g\n",(double)norm);CHKERRQ(ierr);
224         }
225       }
226 
227       if (ipack == 2 && size == 1) {
228         Mat spRHS,spRHST,RHST;
229 
230         ierr = MatTranspose(RHS,MAT_INITIAL_MATRIX,&RHST);CHKERRQ(ierr);
231         ierr = MatConvert(RHST,MATAIJ,MAT_INITIAL_MATRIX,&spRHST);CHKERRQ(ierr);
232         ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr);
233         for (nsolve = 0; nsolve < 2; nsolve++) {
234           ierr = PetscPrintf(PETSC_COMM_WORLD,"   %D-the sparse MatMatSolve \n",nsolve);CHKERRQ(ierr);
235           ierr = MatMatSolve(F,spRHS,X);CHKERRQ(ierr);
236 
237           /* Check the error */
238           ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
239           ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
240           if (norm > tol) {
241             ierr = PetscPrintf(PETSC_COMM_WORLD,"%D-the sparse MatMatSolve: Norm of error %g, nsolve %D\n",nsolve,(double)norm,nsolve);CHKERRQ(ierr);
242           }
243         }
244         ierr = MatDestroy(&spRHST);CHKERRQ(ierr);
245         ierr = MatDestroy(&spRHS);CHKERRQ(ierr);
246         ierr = MatDestroy(&RHST);CHKERRQ(ierr);
247       }
248     }
249 
250     /* Test MatSolve() */
251     if (testMatSolve) {
252       for (nsolve = 0; nsolve < 2; nsolve++) {
253         ierr = VecSetRandom(x,rand);CHKERRQ(ierr);
254         ierr = VecCopy(x,u);CHKERRQ(ierr);
255         ierr = MatMult(A,x,b);CHKERRQ(ierr);
256 
257         ierr = PetscPrintf(PETSC_COMM_WORLD,"   %D-the MatSolve \n",nsolve);CHKERRQ(ierr);
258         ierr = MatSolve(F,b,x);CHKERRQ(ierr);
259 
260         /* Check the error */
261         ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr);  /* u <- (-1.0)x + u */
262         ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
263         if (norm > tol) {
264           PetscReal resi;
265           ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */
266           ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr);  /* u <- (-1.0)b + u */
267           ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr);
268           ierr = PetscPrintf(PETSC_COMM_WORLD,"MatSolve: Norm of error %g, resi %g, numfact %D\n",(double)norm,(double)resi,nfact);CHKERRQ(ierr);
269         }
270       }
271     }
272   }
273 
274   /* Free data structures */
275   ierr = MatDestroy(&A);CHKERRQ(ierr);
276   ierr = MatDestroy(&C);CHKERRQ(ierr);
277   ierr = MatDestroy(&F);CHKERRQ(ierr);
278   ierr = MatDestroy(&X);CHKERRQ(ierr);
279   if (testMatMatSolve) {
280     ierr = MatDestroy(&RHS);CHKERRQ(ierr);
281   }
282 
283   ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
284   ierr = ISDestroy(&perm);CHKERRQ(ierr);
285   ierr = ISDestroy(&iperm);CHKERRQ(ierr);
286   ierr = VecDestroy(&x);CHKERRQ(ierr);
287   ierr = VecDestroy(&b);CHKERRQ(ierr);
288   ierr = VecDestroy(&u);CHKERRQ(ierr);
289   ierr = PetscFinalize();
290   return ierr;
291 }
292 
293 /*TEST
294 
295    test:
296       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
297       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
298       output_file: output/ex125.out
299 
300    test:
301       suffix: mkl_pardiso
302       requires: mkl_pardiso datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
303       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3
304 
305    test:
306       suffix: mumps
307       requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
308       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
309       output_file: output/ex125_mumps_seq.out
310 
311    test:
312       suffix: mumps_2
313       nsize: 3
314       requires: mumps datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
315       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
316       output_file: output/ex125_mumps_par.out
317 
318    test:
319       suffix: superlu_dist
320       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
321       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
322 
323    test:
324       suffix: superlu_dist_2
325       nsize: 3
326       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) superlu_dist
327       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
328       output_file: output/ex125_superlu_dist.out
329 
330    test:
331       suffix: superlu_dist_complex
332       nsize: 3
333       requires: datafilespath superlu_dist complex double !define(PETSC_USE_64BIT_INDICES)
334       args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
335       output_file: output/ex125_superlu_dist_complex.out
336 
337    test:
338       suffix: cusparse
339       requires: cuda datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
340       args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output}
341 
342 TEST*/
343