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,start,end;
15: PetscScalar *array,rval;
16: PetscReal norm,tol=1.e-12;
17: IS perm,iperm;
18: MatFactorInfo info;
19: PetscRandom rand;
20: PetscBool 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 C */
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);
33: /* Assembly C */
34: MatGetOwnershipRange(C,&start,&end);
35: m = end - start;
36: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
37: PetscRandomSetFromOptions(rand);
38: MatGetArray(C,&array);
39: for (i=0; i<m*M; i++){
40: PetscRandomGetValue(rand,&rval);
41: array[i] = rval;
42: }
43: MatRestoreArray(C,&array);
44: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
45: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
46: /* if (!rank) {printf("main, C: \n");}
47: MatView(C,PETSC_VIEWER_STDOUT_WORLD); */
48:
49: /* Create vectors */
50: MatGetLocalSize(C,&m,&n);
51: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix local size m %d must equal n %d",m,n);
52: VecCreate(PETSC_COMM_WORLD,&x);
53: VecSetSizes(x,n,PETSC_DECIDE);
54: VecSetFromOptions(x);
55: VecDuplicate(x,&b);
56: VecDuplicate(x,&u); /* save the true solution */
58: /* Test MatDuplicate() */
59: MatDuplicate(C,MAT_COPY_VALUES,&C1);
60: MatEqual(C,C1,&flg);
61: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Duplicate C1 != C");
63: /* Test LU Factorization */
64: MatGetOrdering(C1,MATORDERINGNATURAL,&perm,&iperm);
65: if (nproc == 1){
66: MatGetFactor(C1,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
67: } else {
68: MatGetFactor(C1,MATSOLVERPLAPACK,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,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
113: } else {
114: MatGetFactor(C,MATSOLVERPLAPACK,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: }