Actual source code: ex125.c
1:
2: static char help[] = "Tests MatSolve and MatMatSolve (interface to superlu_dist).\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,nproc;
15: PetscInt i,m,n,nfact,nsolve,nrhs,k,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, &nproc);
29: /* Determine file from which we read the matrix A */
30: PetscOptionsGetString(PETSC_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) {
40: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
41: }
42:
43: /* Create dense matrix C and X; C holds true solution with identical colums */
44: nrhs = 2;
45: PetscOptionsGetInt(PETSC_NULL,"-nrhs",&nrhs,PETSC_NULL);
46: if (!rank) printf("ex125: nrhs %d\n",nrhs);
47: MatCreate(PETSC_COMM_WORLD,&C);
48: MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
49: MatSetType(C,MATDENSE);
50: MatSetFromOptions(C);
51:
52: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
53: PetscRandomSetFromOptions(rand);
54: MatGetArray(C,&array);
55: for (i=0; i<m; i++){
56: PetscRandomGetValue(rand,&rval);
57: array[i] = rval;
58: }
59: if (nrhs > 1){
60: for (k=1; k<nrhs; k++){
61: for (i=0; i<m; i++){
62: array[m*k+i] = array[i];
63: }
64: }
65: }
66: MatRestoreArray(C,&array);
67: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
69: MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);
70:
71: /* Create vectors */
72: VecCreate(PETSC_COMM_WORLD,&x);
73: VecSetSizes(x,n,PETSC_DECIDE);
74: VecSetFromOptions(x);
75: VecDuplicate(x,&b);
76: VecDuplicate(x,&u); /* save the true solution */
78: /* Test LU Factorization */
79: MatGetOrdering(A,MATORDERINGND,&perm,&iperm);
80: //ISView(perm,PETSC_VIEWER_STDOUT_WORLD);
81: //ISView(perm,PETSC_VIEWER_STDOUT_SELF);
82:
83: PetscOptionsGetInt(PETSC_NULL,"-mat_solver_package",&ipack,PETSC_NULL);
84: switch (ipack){
85: case 0:
86: #ifdef PETSC_HAVE_SUPERLU
87: if (!rank) printf(" SUPERLU LU:\n");
88: MatGetFactor(A,MATSOLVERSUPERLU,MAT_FACTOR_LU,&F);
89: break;
90: #endif
91: case 1:
92: #ifdef PETSC_HAVE_SUPERLU_DIST
93: if (!rank) printf(" SUPERLU_DIST LU:\n");
94: MatGetFactor(A,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&F);
95: break;
96: #endif
97: case 2:
98: #ifdef PETSC_HAVE_MUMPS
99: if (!rank) printf(" MUMPS LU:\n");
100: MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);
101: {
102: /* test mumps options */
103: PetscInt icntl_7 = 5;
104: MatMumpsSetIcntl(F,7,icntl_7);
105: }
106: break;
107: #endif
108: default:
109: if (!rank) printf(" PETSC LU:\n");
110: MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
111: }
113: info.fill = 5.0;
114: MatLUFactorSymbolic(F,A,perm,iperm,&info);
116: for (nfact = 0; nfact < 2; nfact++){
117: if (!rank) printf(" %d-the LU numfactorization \n",nfact);
118: MatLUFactorNumeric(F,A,&info);
120: /* Test MatMatSolve() */
121: if ((ipack == 0 || ipack == 2) && testMatMatSolve){
122: printf(" MaMattSolve() is not implemented for this package. Skip the testing.\n");
123: testMatMatSolve = PETSC_FALSE;
124: }
125: if (testMatMatSolve){
126: if (!nfact){
127: MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
128: } else {
129: MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);
130: }
131: for (nsolve = 0; nsolve < 2; nsolve++){
132: if (!rank) printf(" %d-the MatMatSolve \n",nsolve);
133: MatMatSolve(F,RHS,X);
134:
135: /* Check the error */
136: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
137: MatNorm(X,NORM_FROBENIUS,&norm);
138: if (norm > tol){
139: if (!rank){
140: PetscPrintf(PETSC_COMM_SELF,"1st MatMatSolve: Norm of error %g, nsolve %d\n",norm,nsolve);
141: }
142: }
143: }
144: }
146: /* Test MatSolve() */
147: if (testMatSolve){
148: for (nsolve = 0; nsolve < 2; nsolve++){
149: VecGetArray(x,&array);
150: for (i=0; i<m; i++){
151: PetscRandomGetValue(rand,&rval);
152: array[i] = rval;
153: }
154: VecRestoreArray(x,&array);
155: VecCopy(x,u);
156: MatMult(A,x,b);
158: if (!rank) printf(" %d-the MatSolve \n",nsolve);
159: MatSolve(F,b,x);
160:
161: /* Check the error */
162: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
163: VecNorm(u,NORM_2,&norm);
164: if (norm > tol){
165: MatMult(A,x,u); /* u = A*x */
166: PetscReal resi;
167: VecAXPY(u,-1.0,b); /* u <- (-1.0)b + u */
168: VecNorm(u,NORM_2,&resi);
169: if (!rank){
170: PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g, resi %g, LU numfact %d\n",norm,resi,nfact);
171: }
172: }
173: }
174: }
176: /* Test MatMatSolve() */
177: if (testMatMatSolve){
178: for (nsolve = 0; nsolve < 2; nsolve++){
179: if (!rank) printf(" %d-the MatMatSolve \n",nsolve);
180: MatMatSolve(F,RHS,X);
181:
182: /* Check the error */
183: MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
184: MatNorm(X,NORM_FROBENIUS,&norm);
185: if (norm > tol){
186: if (!rank){
187: PetscPrintf(PETSC_COMM_SELF,"2nd MatMatSolve: Norm of error %g, nsolve %d\n",norm,nsolve);
188: }
189: }
190: }
191: }
192: }
193:
194: /* Free data structures */
195: MatDestroy(&A);
196: MatDestroy(&C);
197: MatDestroy(&F);
198: MatDestroy(&X);
199: if (testMatMatSolve){
200: MatDestroy(&RHS);
201: }
202:
203: PetscRandomDestroy(&rand);
204: ISDestroy(&perm);
205: ISDestroy(&iperm);
206: VecDestroy(&x);
207: VecDestroy(&b);
208: VecDestroy(&u);
209: PetscFinalize();
210: return 0;
211: }