Actual source code: ex103.c
petsc-3.3-p6 2013-02-11
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);
32: MatSetUp(C);
34: /* Assembly C */
35: MatGetOwnershipRange(C,&start,&end);
36: m = end - start;
37: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
38: PetscRandomSetFromOptions(rand);
39: MatGetArray(C,&array);
40: for (i=0; i<m*M; i++){
41: PetscRandomGetValue(rand,&rval);
42: array[i] = rval;
43: }
44: MatRestoreArray(C,&array);
45: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
46: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
47: /* if (!rank) {printf("main, C: \n");}
48: MatView(C,PETSC_VIEWER_STDOUT_WORLD); */
49:
50: /* Create vectors */
51: MatGetLocalSize(C,&m,&n);
52: if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix local size m %d must equal n %d",m,n);
53: VecCreate(PETSC_COMM_WORLD,&x);
54: VecSetSizes(x,n,PETSC_DECIDE);
55: VecSetFromOptions(x);
56: VecDuplicate(x,&b);
57: VecDuplicate(x,&u); /* save the true solution */
59: /* Test MatDuplicate() */
60: MatDuplicate(C,MAT_COPY_VALUES,&C1);
61: MatEqual(C,C1,&flg);
62: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Duplicate C1 != C");
64: /* Test LU Factorization */
65: MatGetOrdering(C1,MATORDERINGNATURAL,&perm,&iperm);
66: if (nproc == 1){
67: MatGetFactor(C1,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
68: } else {
69: MatGetFactor(C1,MATSOLVERPLAPACK,MAT_FACTOR_LU,&F);
70: }
71: MatLUFactorSymbolic(F,C1,perm,iperm,&info);
73: for (nfact = 0; nfact < 2; nfact++){
74: if (!rank) printf(" LU nfact %d\n",nfact);
75: MatLUFactorNumeric(F,C1,&info);
77: /* Test MatSolve() */
78: for (nsolve = 0; nsolve < 5; nsolve++){
79: VecGetArray(x,&array);
80: for (i=0; i<m; i++){
81: PetscRandomGetValue(rand,&rval);
82: array[i] = rval;
83: }
84: VecRestoreArray(x,&array);
85: VecCopy(x,u);
86: MatMult(C,x,b);
88: MatSolve(F,b,x);
90: /* Check the error */
91: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
92: VecNorm(u,NORM_2,&norm);
93: if (norm > tol){
94: if (!rank){
95: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, LU nfact %d\n",norm,nfact);
96: }
97: }
98: }
99: }
100: MatDestroy(&C1);
101: MatDestroy(&F);
103: /* Test Cholesky Factorization */
104: MatTranspose(C,MAT_INITIAL_MATRIX,&C1); /* C1 = C^T */
105: MatAXPY(C,1.0,C1,SAME_NONZERO_PATTERN); /* make C symmetric: C <- C + C^T */
106: MatShift(C,M); /* make C positive definite */
107: MatDestroy(&C1);
108:
109: MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
110: MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
111:
112: if (nproc == 1){
113: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
114: } else {
115: MatGetFactor(C,MATSOLVERPLAPACK,MAT_FACTOR_CHOLESKY,&F);
116: }
117: MatCholeskyFactorSymbolic(F,C,perm,&info);
118: for (nfact = 0; nfact < 2; nfact++){
119: if (!rank) printf(" Cholesky nfact %d\n",nfact);
120: MatCholeskyFactorNumeric(F,C,&info);
122: /* Test MatSolve() */
123: for (nsolve = 0; nsolve < 5; nsolve++){
124: VecGetArray(x,&array);
125: for (i=0; i<m; i++){
126: PetscRandomGetValue(rand,&rval);
127: array[i] = rval;
128: }
129: VecRestoreArray(x,&array);
130: VecCopy(x,u);
131: MatMult(C,x,b);
133: MatSolve(F,b,x);
135: /* Check the error */
136: VecAXPY(u,-1.0,x); /* u <- (-1.0)x + u */
137: VecNorm(u,NORM_2,&norm);
138: if (norm > tol){
139: if (!rank){
140: PetscPrintf(PETSC_COMM_SELF,"Error: Norm of error %g, Cholesky nfact %d\n",norm,nfact);
141: }
142: }
143: }
144: }
145: MatDestroy(&F);
147: /* Free data structures */
148: PetscRandomDestroy(&rand);
149: ISDestroy(&perm);
150: ISDestroy(&iperm);
151: VecDestroy(&x);
152: VecDestroy(&b);
153: VecDestroy(&u);
154: MatDestroy(&C);
156: PetscFinalize();
157: return 0;
158: }