Actual source code: ex1.c
2: static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a sequential dense matrix. \n\
3: For MATSEQDENSE matrix, the factorization is just a thin wrapper to LAPACK \n\n";
5: #include <petscmat.h>
9: int main(int argc,char **argv)
10: {
11: Mat mat,F,RHS,SOLU;
12: MatInfo info;
14: PetscInt m = 10,n = 10,i,j,rstart,rend,nrhs=2;
15: PetscScalar value = 1.0;
16: Vec x,y,b,ytmp;
17: PetscReal norm;
18: PetscMPIInt size;
19: PetscScalar *rhs_array,*solu_array;
20: PetscRandom rand;
21: PetscScalar *array,rval;
23: PetscInitialize(&argc,&argv,(char*) 0,help);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
27: /* create single vectors */
28: VecCreate(PETSC_COMM_WORLD,&y);
29: VecSetSizes(y,PETSC_DECIDE,m);
30: VecSetFromOptions(y);
31: VecDuplicate(y,&x);
32: VecDuplicate(y,&ytmp);
33: VecSet(x,value);
34: VecCreate(PETSC_COMM_WORLD,&b);
35: VecSetSizes(b,PETSC_DECIDE,n);
36: VecSetFromOptions(b);
38: /* create multiple vectors RHS and SOLU */
39: MatCreate(PETSC_COMM_WORLD,&RHS);
40: MatSetSizes(RHS,PETSC_DECIDE,PETSC_DECIDE,n,nrhs);
41: MatSetType(RHS,MATDENSE);
42: MatSetFromOptions(RHS);
43:
44: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
45: PetscRandomSetFromOptions(rand);
46: MatGetArray(RHS,&array);
47: for (j=0; j<nrhs; j++){
48: for (i=0; i<n; i++){
49: PetscRandomGetValue(rand,&rval);
50: array[n*j+i] = rval;
51: }
52: }
53: MatRestoreArray(RHS,&array);
54:
55: MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&SOLU);
57: /* create matrix */
58: MatCreateSeqDense(PETSC_COMM_WORLD,m,n,PETSC_NULL,&mat);
59: MatGetOwnershipRange(mat,&rstart,&rend);
60: for (i=rstart; i<rend; i++) {
61: value = (PetscReal)i+1;
62: MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
63: }
64: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
67: MatGetInfo(mat,MAT_LOCAL,&info);
68: PetscPrintf(PETSC_COMM_WORLD,"matrix nonzeros = %D, allocated nonzeros = %D\n",
69: (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
71: /* Cholesky factorization - perm and factinfo are ignored by LAPACK */
72: /* in-place Cholesky */
73: MatMult(mat,x,b);
74: MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&F);
75: MatCholeskyFactor(F,0,0);
76: MatSolve(F,b,y);
77: MatDestroy(&F);
78: value = -1.0; VecAXPY(y,value,x);
79: VecNorm(y,NORM_2,&norm);
80: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
81: /* out-place Cholesky */
82: MatGetFactor(mat,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&F);
83: MatCholeskyFactorSymbolic(F,mat,0,0);
84: MatCholeskyFactorNumeric(F,mat,0);
85: MatSolve(F,b,y);
86: value = -1.0; VecAXPY(y,value,x);
87: VecNorm(y,NORM_2,&norm);
88: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for Cholesky %A\n",norm);
89: MatDestroy(&F);
91: /* LU factorization - perms and factinfo are ignored by LAPACK */
92: i = m-1; value = 1.0;
93: MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
94: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
96: MatMult(mat,x,b);
97: MatConvert(mat,MATSAME,MAT_INITIAL_MATRIX,&F);
99: /* in-place LU */
100: MatLUFactor(F,0,0,0);
101: MatSolve(F,b,y);
102: value = -1.0; VecAXPY(y,value,x);
103: VecNorm(y,NORM_2,&norm);
104: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
106: MatMatSolve(F,RHS,SOLU);
107: for (j=0; j<nrhs; j++){
108: MatGetArray(SOLU,&solu_array);
109: MatGetArray(RHS,&rhs_array);
110: VecPlaceArray(y,solu_array);
111: VecPlaceArray(b,rhs_array);
113: MatMult(mat,y,ytmp);
114: VecAXPY(ytmp,-1.0,b); /* ytmp = mat*SOLU[:,j] - RHS[:,j] */
115: VecNorm(ytmp,NORM_2,&norm);
116: if (norm > 1.e-12){
117: PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for LU %A\n",norm);
118: }
119:
120: VecResetArray(b);
121: VecResetArray(y);
122: MatRestoreArray(RHS,&rhs_array);
123: MatRestoreArray(SOLU,&solu_array);
124: }
126: MatDestroy(&F);
127: /* out-place LU */
128: MatGetFactor(mat,MATSOLVERPETSC,MAT_FACTOR_LU,&F);
129: MatLUFactorSymbolic(F,mat,0,0,0);
130: MatLUFactorNumeric(F,mat,0);
131: MatSolve(F,b,y);
132: value = -1.0; VecAXPY(y,value,x);
133: VecNorm(y,NORM_2,&norm);
134: PetscPrintf(PETSC_COMM_WORLD,"Norm of error for LU %A\n",norm);
136: /* free space */
137: MatDestroy(&F);
138: MatDestroy(&mat);
139: MatDestroy(&RHS);
140: MatDestroy(&SOLU);
141: PetscRandomDestroy(&rand);
142: VecDestroy(&x);
143: VecDestroy(&b);
144: VecDestroy(&y);
145: VecDestroy(&ytmp);
146: PetscFinalize();
147: return 0;
148: }
149: