Actual source code: ex1.c
petsc-3.14.6 2021-03-30
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\
4: For MATSEQDENSECUDA, it uses cusolverDn routines \n\n";
6: #include <petscmat.h>
8: int main(int argc,char **argv)
9: {
10: Mat mat,F,RHS,SOLU;
11: MatInfo info;
13: PetscInt n = 10,i,j,rstart,rend,nrhs=2;
14: PetscScalar value = 1.0;
15: Vec x,y,b,ytmp;
16: IS perm;
17: PetscReal norm,tol=PETSC_SMALL;
18: PetscMPIInt size;
19: PetscRandom rand;
20: char solver[64];
21: PetscBool inplace,full = PETSC_FALSE, ldl = PETSC_TRUE;
23: PetscInitialize(&argc,&argv,(char*) 0,help);if (ierr) return ierr;
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!");
26: PetscStrcpy(solver,"petsc");
27: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
29: PetscOptionsGetBool(NULL,NULL,"-ldl",&ldl,NULL);
30: PetscOptionsGetBool(NULL,NULL,"-full",&full,NULL);
31: PetscOptionsGetString(NULL,NULL,"-solver_type",solver,sizeof(solver),NULL);
33: /* create multiple vectors RHS and SOLU */
34: MatCreate(PETSC_COMM_WORLD,&RHS);
35: MatSetSizes(RHS,PETSC_DECIDE,PETSC_DECIDE,n,nrhs);
36: MatSetType(RHS,MATDENSE);
37: MatSetOptionsPrefix(RHS,"rhs_");
38: MatSetFromOptions(RHS);
39: MatSeqDenseSetPreallocation(RHS,NULL);
41: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
42: PetscRandomSetFromOptions(rand);
43: MatSetRandom(RHS,rand);
45: MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&SOLU);
47: /* create matrix */
48: MatCreate(PETSC_COMM_WORLD,&mat);
49: MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,n,n);
50: MatSetType(mat,MATDENSE);
51: MatSetFromOptions(mat);
52: MatSetUp(mat);
53: MatGetOwnershipRange(mat,&rstart,&rend);
54: if (!full) {
55: for (i=rstart; i<rend; i++) {
56: value = (PetscReal)i+1;
57: MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);
58: }
59: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
60: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
61: } else {
62: Mat T;
64: MatSetRandom(mat,rand);
65: MatMatTransposeMult(mat,mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);
66: MatDestroy(&mat);
67: mat = T;
68: }
70: /* create single vectors */
71: MatCreateVecs(mat,&x,&b);
72: VecDuplicate(x,&y);
73: VecDuplicate(y,&ytmp);
74: VecSet(x,value);
76: /* Only SeqDense* support in-place factorizations and NULL permutations */
77: PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQDENSE,&inplace);
78: MatGetLocalSize(mat,&i,NULL);
79: MatGetOwnershipRange(mat,&j,NULL);
80: ISCreateStride(PETSC_COMM_WORLD,i,j,1,&perm);
82: MatGetInfo(mat,MAT_LOCAL,&info);
83: PetscPrintf(PETSC_COMM_WORLD,"matrix nonzeros = %D, allocated nonzeros = %D\n",
84: (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);
85: MatMult(mat,x,b);
87: /* Cholesky factorization - perm and factinfo are ignored by LAPACK */
88: /* in-place Cholesky */
89: if (inplace) {
90: Mat RHS2;
92: MatDuplicate(mat,MAT_COPY_VALUES,&F);
93: if (!ldl) { MatSetOption(F,MAT_SPD,PETSC_TRUE); }
94: MatCholeskyFactor(F,perm,0);
95: MatSolve(F,b,y);
96: VecAXPY(y,-1.0,x);
97: VecNorm(y,NORM_2,&norm);
98: if (norm > tol) {
99: PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place Cholesky %g\n",(double)norm);
100: }
102: MatMatSolve(F,RHS,SOLU);
103: MatMatMult(mat,SOLU,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RHS2);
104: MatAXPY(RHS,-1.0,RHS2,SAME_NONZERO_PATTERN);
105: MatNorm(RHS,NORM_FROBENIUS,&norm);
106: if (norm > tol) {
107: PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for in-place Cholesky (MatMatSolve) %g\n",(double)norm);
108: }
109: MatDestroy(&F);
110: MatDestroy(&RHS2);
111: }
113: /* out-of-place Cholesky */
114: MatGetFactor(mat,solver,MAT_FACTOR_CHOLESKY,&F);
115: if (!ldl) { MatSetOption(F,MAT_SPD,PETSC_TRUE); }
116: MatCholeskyFactorSymbolic(F,mat,perm,0);
117: MatCholeskyFactorNumeric(F,mat,0);
118: MatSolve(F,b,y);
119: VecAXPY(y,-1.0,x);
120: VecNorm(y,NORM_2,&norm);
121: if (norm > tol) {
122: PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place Cholesky %g\n",(double)norm);
123: }
124: MatDestroy(&F);
126: /* LU factorization - perms and factinfo are ignored by LAPACK */
127: i = n-1;
128: MatZeroRows(mat,1,&i,-1.0,NULL,NULL);
129: MatMult(mat,x,b);
131: /* in-place LU */
132: if (inplace) {
133: Mat RHS2;
135: MatDuplicate(mat,MAT_COPY_VALUES,&F);
136: MatLUFactor(F,perm,perm,0);
137: MatSolve(F,b,y);
138: VecAXPY(y,-1.0,x);
139: VecNorm(y,NORM_2,&norm);
140: if (norm > tol) {
141: PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place LU %g\n",(double)norm);
142: }
143: MatMatSolve(F,RHS,SOLU);
144: MatMatMult(mat,SOLU,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RHS2);
145: MatAXPY(RHS,-1.0,RHS2,SAME_NONZERO_PATTERN);
146: MatNorm(RHS,NORM_FROBENIUS,&norm);
147: if (norm > tol) {
148: PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for in-place LU (MatMatSolve) %g\n",(double)norm);
149: }
150: MatDestroy(&F);
151: MatDestroy(&RHS2);
152: }
154: /* out-of-place LU */
155: MatGetFactor(mat,solver,MAT_FACTOR_LU,&F);
156: MatLUFactorSymbolic(F,mat,perm,perm,0);
157: MatLUFactorNumeric(F,mat,0);
158: MatSolve(F,b,y);
159: VecAXPY(y,-1.0,x);
160: VecNorm(y,NORM_2,&norm);
161: if (norm > tol) {
162: PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place LU %g\n",(double)norm);
163: }
165: /* free space */
166: ISDestroy(&perm);
167: MatDestroy(&F);
168: MatDestroy(&mat);
169: MatDestroy(&RHS);
170: MatDestroy(&SOLU);
171: PetscRandomDestroy(&rand);
172: VecDestroy(&x);
173: VecDestroy(&b);
174: VecDestroy(&y);
175: VecDestroy(&ytmp);
176: PetscFinalize();
177: return ierr;
178: }
182: /*TEST
184: test:
186: test:
187: requires: cuda
188: suffix: seqdensecuda
189: args: -mat_type seqdensecuda -rhs_mat_type seqdensecuda -ldl 0 -solver_type {{petsc cuda}}
190: output_file: output/ex1_1.out
192: test:
193: requires: cuda
194: suffix: seqdensecuda_seqaijcusparse
195: args: -mat_type seqaijcusparse -rhs_mat_type seqdensecuda
196: output_file: output/ex1_2.out
198: test:
199: requires: cuda viennacl
200: suffix: seqdensecuda_seqaijviennacl
201: args: -mat_type seqaijviennacl -rhs_mat_type seqdensecuda
202: output_file: output/ex1_2.out
204: TEST*/