Actual source code: ex27.c
petsc-3.3-p6 2013-02-11
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: Test MatMatSolve(). Input parameters include\n\
4: -f <input_file> : file to load \n\n";
6: /*
7: Usage:
8: ex27 -f0 <mat_binaryfile>
9: */
11: #include <petscksp.h>
12: extern PetscErrorCode PCShellApply_Matinv(PC,Vec,Vec);
16: int main(int argc,char **args)
17: {
18: KSP ksp;
19: Mat A,B,F,X;
20: Vec x,b,u; /* approx solution, RHS, exact solution */
21: PetscViewer fd; /* viewer */
22: char file[1][PETSC_MAX_PATH_LEN]; /* input file name */
23: PetscBool flg;
25: PetscInt M,N,i,its;
26: PetscReal norm;
27: PetscScalar val=1.0;
28: PetscMPIInt size;
29: PC pc;
31: PetscInitialize(&argc,&args,(char *)0,help);
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
33: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
35: /* Read matrix and right-hand-side vector */
36: PetscOptionsGetString(PETSC_NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
37: if (!flg) {
38: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
39: }
41: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);
42: MatCreate(PETSC_COMM_WORLD,&A);
43: MatSetType(A,MATAIJ);
44: MatLoad(A,fd);
45: VecCreate(PETSC_COMM_WORLD,&b);
46: VecLoad(b,fd);
47: PetscViewerDestroy(&fd);
49: /*
50: If the loaded matrix is larger than the vector (due to being padded
51: to match the block size of the system), then create a new padded vector.
52: */
53: {
54: PetscInt m,n,j,mvec,start,end,indx;
55: Vec tmp;
56: PetscScalar *bold;
58: /* Create a new vector b by padding the old one */
59: MatGetLocalSize(A,&m,&n);
60: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
61: VecCreate(PETSC_COMM_WORLD,&tmp);
62: VecSetSizes(tmp,m,PETSC_DECIDE);
63: VecSetFromOptions(tmp);
64: VecGetOwnershipRange(b,&start,&end);
65: VecGetLocalSize(b,&mvec);
66: VecGetArray(b,&bold);
67: for (j=0; j<mvec; j++) {
68: indx = start+j;
69: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
70: }
71: VecRestoreArray(b,&bold);
72: VecDestroy(&b);
73: VecAssemblyBegin(tmp);
74: VecAssemblyEnd(tmp);
75: b = tmp;
76: }
77: VecDuplicate(b,&x);
78: VecDuplicate(b,&u);
79: VecSet(x,0.0);
81: /* Create dense matric B and X. Set B as an identity matrix */
82: MatGetSize(A,&M,&N);
83: MatCreate(MPI_COMM_SELF,&B);
84: MatSetSizes(B,M,N,M,N);
85: MatSetType(B,MATSEQDENSE);
86: MatSeqDenseSetPreallocation(B,PETSC_NULL);
87: for (i=0; i<M; i++){
88: MatSetValues(B,1,&i,1,&i,&val,INSERT_VALUES);
89: }
90: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
93: MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X);
94:
95: /* Compute X=inv(A) by MatMatSolve() */
96: KSPCreate(PETSC_COMM_WORLD,&ksp);
97: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
98: KSPGetPC(ksp,&pc);
99: PCSetType(pc,PCLU);
100: KSPSetUp(ksp);
101: PCFactorGetMatrix(pc,&F);
102: MatMatSolve(F,B,X);
103: MatDestroy(&B);
105: /* Now, set X=inv(A) as a preconditioner */
106: PCSetType(pc,PCSHELL);
107: PCShellSetContext(pc,(void *)X);
108: PCShellSetApply(pc,PCShellApply_Matinv);
109: KSPSetFromOptions(ksp);
111: /* Sove preconditioned system A*x = b */
112: KSPSolve(ksp,b,x);
113: KSPGetIterationNumber(ksp,&its);
115: /* Check error */
116: MatMult(A,x,u);
117: VecAXPY(u,-1.0,b);
118: VecNorm(u,NORM_2,&norm);
119: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
120: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
121:
122: /* Free work space. */
123: MatDestroy(&X);
124: MatDestroy(&A); VecDestroy(&b);
125: VecDestroy(&u); VecDestroy(&x);
126: KSPDestroy(&ksp);
127: PetscFinalize();
128: return 0;
129: }
131: PetscErrorCode PCShellApply_Matinv(PC pc,Vec xin,Vec xout)
132: {
134: Mat X;
137: PCShellGetContext(pc,(void**)&X);
138: MatMult(X,xin,xout);
139: return(0);
140: }