Actual source code: ex103.c
1:
2: static char help[] = "Tests PLAPACK interface.\n\n";
4: #include petscmat.h
8: int main(int argc,char **args)
9: {
10: Mat C,C1,F;
11: Vec u,x,b;
13: PetscMPIInt rank,nproc;
14: PetscInt i,M = 10,m,n,nfact,nsolve;
15: PetscScalar *array,rval;
16: PetscReal norm,tol=1.e-12;
17: IS perm,iperm;
18: MatFactorInfo info;
19: PetscRandom rand;
20: PetscTruth flg;
22: PetscInitialize(&argc,&args,(char *)0,help);
23: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
24: MPI_Comm_size(PETSC_COMM_WORLD, &nproc);
26: /* Create matrix and vectors */
27: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
28: MatCreate(PETSC_COMM_WORLD,&C);
29: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,M,M);
30: MatSetType(C,MATDENSE);
31: MatSetFromOptions(C);
32:
33: MatGetLocalSize(C,&m,&n);
34: if (m != n) SETERRQ2(PETSC_ERR_ARG_WRONG,"Matrix local size m %d must equal n %d",m,n);
36: VecCreate(PETSC_COMM_WORLD,&x);
37: VecSetSizes(x,n,PETSC_DECIDE);
38: VecSetFromOptions(x);
39: VecDuplicate(x,&b);
40: VecDuplicate(x,&u); /* save the true solution */
42: /* Assembly */
43: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
44: PetscRandomSetFromOptions(rand);
45: MatGetArray(C,&array);
46: for (i=0; i<m*M; i++){
47: PetscRandomGetValue(rand,&rval);
48: array[i] = rval;
49: }
50: MatRestoreArray(C,&array);
51: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
52: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
53: /*if (!rank) {printf("main, C: \n");}
54: MatView(C,PETSC_VIEWER_STDOUT_WORLD); */
56: /* Test MatDuplicate() */
57: MatDuplicate(C,MAT_COPY_VALUES,&C1);
58: MatEqual(C,C1,&flg);
59: if (!flg){
60: SETERRQ(PETSC_ERR_ARG_WRONG,"Duplicate C1 != C");
61: }
63: /* Test LU Factorization */
64: MatGetOrdering(C1,MATORDERING_NATURAL,&perm,&iperm);
65: if (nproc == 1){
66: MatGetFactor(C1,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&F);
67: } else {
68: MatGetFactor(C1,MAT_SOLVER_PLAPACK,MAT_FACTOR_LU,&F);
69: }
70: MatLUFactorSymbolic(F,C1,perm,iperm,&info);
72: for (nfact = 0; nfact < 2; nfact++){
73: if (!rank) printf(" LU nfact %d\n",nfact);
74: MatLUFactorNumeric(F,C1,&info);
76: /* Test MatSolve() */
77: for (nsolve = 0; nsolve < 5; nsolve++){
78: VecGetArray(x,&array);
79: for (i=0; i<m; i++){
80: PetscRandomGetValue(rand,&rval);
81: array[i] = rval;
82: }
83: VecRestoreArray(x,&array);
84: VecCopy(x,u);
85: MatMult(C,x,b);
87: MatSolve(F,b,x);
89: /* Check the error */
90: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
91: VecNorm(u,NORM_2,&norm);
92: if (norm > tol){
93: if (!rank){
94: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, LU nfact %d\n",norm,nfact);
95: }
96: }
97: }
98: }
99: MatDestroy(C1);
100: MatDestroy(F);
102: /* Test Cholesky Factorization */
103: MatTranspose(C,MAT_INITIAL_MATRIX,&C1); /* C1 = C^T */
104: MatAXPY(C,1.0,C1,SAME_NONZERO_PATTERN); /* make C symmetric: C <- C + C^T */
105: MatShift(C,M); /* make C positive definite */
106: MatDestroy(C1);
107:
108: MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
109: MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
110:
111: if (nproc == 1){
112: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&F);
113: } else {
114: MatGetFactor(C,MAT_SOLVER_PLAPACK,MAT_FACTOR_CHOLESKY,&F);
115: }
116: MatCholeskyFactorSymbolic(F,C,perm,&info);
117: for (nfact = 0; nfact < 2; nfact++){
118: if (!rank) printf(" Cholesky nfact %d\n",nfact);
119: MatCholeskyFactorNumeric(F,C,&info);
121: /* Test MatSolve() */
122: for (nsolve = 0; nsolve < 5; nsolve++){
123: VecGetArray(x,&array);
124: for (i=0; i<m; i++){
125: PetscRandomGetValue(rand,&rval);
126: array[i] = rval;
127: }
128: VecRestoreArray(x,&array);
129: VecCopy(x,u);
130: MatMult(C,x,b);
132: MatSolve(F,b,x);
134: /* Check the error */
135: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
136: VecNorm(u,NORM_2,&norm);
137: if (norm > tol){
138: if (!rank){
139: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, Cholesky nfact %d\n",norm,nfact);
140: }
141: }
142: }
143: }
144: MatDestroy(F);
146: /* Free data structures */
147: PetscRandomDestroy(rand);
148: ISDestroy(perm);
149: ISDestroy(iperm);
150: VecDestroy(x);
151: VecDestroy(b);
152: VecDestroy(u);
153: MatDestroy(C);
155: PetscFinalize();
156: return 0;
157: }