xref: /petsc/src/mat/tests/ex1.c (revision 4905a7bc61a644ac28a555b575668251734ce1fa) 
1 
2 static char help[] = "Tests LU, Cholesky, and QR 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";
5 
6 #include <petscmat.h>
7 
8 static PetscErrorCode createMatsAndVecs(PetscInt m, PetscInt n, PetscInt nrhs, PetscBool full, Mat *_mat, Mat *_RHS, Mat *_SOLU, Vec *_x, Vec *_y, Vec *_b)
9 {
10   PetscRandom    rand;
11   Mat            mat,RHS,SOLU;
12   PetscInt       rstart, rend;
13   PetscInt       cstart, cend;
14   PetscScalar    value = 1.0;
15   Vec            x, y, b;
16   PetscErrorCode ierr;
17 
18   PetscFunctionBegin;
19   /* create multiple vectors RHS and SOLU */
20   ierr = MatCreate(PETSC_COMM_WORLD,&RHS);CHKERRQ(ierr);
21   ierr = MatSetSizes(RHS,PETSC_DECIDE,PETSC_DECIDE,m,nrhs);CHKERRQ(ierr);
22   ierr = MatSetType(RHS,MATDENSE);CHKERRQ(ierr);
23   ierr = MatSetOptionsPrefix(RHS,"rhs_");CHKERRQ(ierr);
24   ierr = MatSetFromOptions(RHS);CHKERRQ(ierr);
25   ierr = MatSeqDenseSetPreallocation(RHS,NULL);CHKERRQ(ierr);
26 
27   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr);
28   ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr);
29   ierr = MatSetRandom(RHS,rand);CHKERRQ(ierr);
30 
31   if (m == n) {
32     ierr = MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&SOLU);CHKERRQ(ierr);
33   } else {
34     ierr = MatCreate(PETSC_COMM_WORLD,&SOLU);CHKERRQ(ierr);
35     ierr = MatSetSizes(SOLU,PETSC_DECIDE,PETSC_DECIDE,n,nrhs);CHKERRQ(ierr);
36     ierr = MatSetType(SOLU,MATDENSE);CHKERRQ(ierr);
37     ierr = MatSeqDenseSetPreallocation(SOLU,NULL);CHKERRQ(ierr);
38   }
39   ierr = MatSetRandom(SOLU,rand);CHKERRQ(ierr);
40 
41   /* create matrix */
42   ierr = MatCreate(PETSC_COMM_WORLD,&mat);CHKERRQ(ierr);
43   ierr = MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
44   ierr = MatSetType(mat,MATDENSE);CHKERRQ(ierr);
45   ierr = MatSetFromOptions(mat);CHKERRQ(ierr);
46   ierr = MatSetUp(mat);CHKERRQ(ierr);
47   ierr = MatGetOwnershipRange(mat,&rstart,&rend);CHKERRQ(ierr);
48   ierr = MatGetOwnershipRangeColumn(mat,&cstart,&cend);CHKERRQ(ierr);
49   if (!full) {
50     for (PetscInt i=rstart; i<rend; i++) {
51       if (m == n) {
52         value = (PetscReal)i+1;
53         ierr  = MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr);
54       } else {
55         for (PetscInt j = cstart; j < cend; j++) {
56           value = ((PetscScalar)i+1.)/(PetscSqr(i - j) + 1.);
57           ierr  = MatSetValues(mat,1,&i,1,&j,&value,INSERT_VALUES);CHKERRQ(ierr);
58         }
59       }
60     }
61     ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
62     ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
63   } else {
64     ierr = MatSetRandom(mat,rand);CHKERRQ(ierr);
65     if (m == n) {
66       Mat T;
67 
68       ierr = MatMatTransposeMult(mat,mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr);
69       ierr = MatDestroy(&mat);CHKERRQ(ierr);
70       mat  = T;
71     }
72   }
73 
74   /* create single vectors */
75   ierr  = MatCreateVecs(mat,&x,&b);CHKERRQ(ierr);
76   ierr  = VecDuplicate(x,&y);CHKERRQ(ierr);
77   ierr  = VecSet(x,value);CHKERRQ(ierr);
78   ierr  = PetscRandomDestroy(&rand);CHKERRQ(ierr);
79   *_mat  = mat;
80   *_RHS  = RHS;
81   *_SOLU = SOLU;
82   *_x    = x;
83   *_y    = y;
84   *_b    = b;
85   PetscFunctionReturn(0);
86 }
87 
88 int main(int argc,char **argv)
89 {
90   Mat            mat,F,RHS,SOLU;
91   MatInfo        info;
92   PetscErrorCode ierr;
93   PetscInt       m = 15, n = 10,i,j,nrhs=2;
94   Vec            x,y,b,ytmp;
95   IS             perm;
96   PetscReal      norm,tol=PETSC_SMALL;
97   PetscMPIInt    size;
98   char           solver[64];
99   PetscBool      inplace,full = PETSC_FALSE, ldl = PETSC_TRUE;
100 
101   ierr = PetscInitialize(&argc,&argv,(char*) 0,help);if (ierr) return ierr;
102   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
103   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
104   ierr = PetscStrcpy(solver,"petsc");CHKERRQ(ierr);
105   ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr);
106   ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
107   ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr);
108   ierr = PetscOptionsGetBool(NULL,NULL,"-ldl",&ldl,NULL);CHKERRQ(ierr);
109   ierr = PetscOptionsGetBool(NULL,NULL,"-full",&full,NULL);CHKERRQ(ierr);
110   ierr = PetscOptionsGetString(NULL,NULL,"-solver_type",solver,sizeof(solver),NULL);CHKERRQ(ierr);
111 
112   ierr = createMatsAndVecs(n, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b);CHKERRQ(ierr);
113   ierr = VecDuplicate(y,&ytmp);CHKERRQ(ierr);
114 
115   /* Only SeqDense* support in-place factorizations and NULL permutations */
116   ierr = PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQDENSE,&inplace);CHKERRQ(ierr);
117   ierr = MatGetLocalSize(mat,&i,NULL);CHKERRQ(ierr);
118   ierr = MatGetOwnershipRange(mat,&j,NULL);CHKERRQ(ierr);
119   ierr = ISCreateStride(PETSC_COMM_WORLD,i,j,1,&perm);CHKERRQ(ierr);
120 
121   ierr = MatGetInfo(mat,MAT_LOCAL,&info);CHKERRQ(ierr);
122   ierr = PetscPrintf(PETSC_COMM_WORLD,"matrix nonzeros = %D, allocated nonzeros = %D\n",
123                      (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);CHKERRQ(ierr);
124   ierr = MatMult(mat,x,b);CHKERRQ(ierr);
125 
126   /* Cholesky factorization - perm and factinfo are ignored by LAPACK */
127   /* in-place Cholesky */
128   if (inplace) {
129     Mat RHS2;
130 
131     ierr = MatDuplicate(mat,MAT_COPY_VALUES,&F);CHKERRQ(ierr);
132     if (!ldl) { ierr = MatSetOption(F,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); }
133     ierr = MatCholeskyFactor(F,perm,0);CHKERRQ(ierr);
134     ierr = MatSolve(F,b,y);CHKERRQ(ierr);
135     ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
136     ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
137     if (norm > tol) {
138       ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place Cholesky %g\n",(double)norm);CHKERRQ(ierr);
139     }
140 
141     ierr = MatMatSolve(F,RHS,SOLU);CHKERRQ(ierr);
142     ierr = MatMatMult(mat,SOLU,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RHS2);CHKERRQ(ierr);
143     ierr = MatAXPY(RHS,-1.0,RHS2,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
144     ierr = MatNorm(RHS,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
145     if (norm > tol) {
146       ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for in-place Cholesky (MatMatSolve) %g\n",(double)norm);CHKERRQ(ierr);
147     }
148     ierr = MatDestroy(&F);CHKERRQ(ierr);
149     ierr = MatDestroy(&RHS2);CHKERRQ(ierr);
150   }
151 
152   /* out-of-place Cholesky */
153   ierr = MatGetFactor(mat,solver,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr);
154   if (!ldl) { ierr = MatSetOption(F,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); }
155   ierr = MatCholeskyFactorSymbolic(F,mat,perm,0);CHKERRQ(ierr);
156   ierr = MatCholeskyFactorNumeric(F,mat,0);CHKERRQ(ierr);
157   ierr = MatSolve(F,b,y);CHKERRQ(ierr);
158   ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
159   ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
160   if (norm > tol) {
161     ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place Cholesky %g\n",(double)norm);CHKERRQ(ierr);
162   }
163   ierr = MatDestroy(&F);CHKERRQ(ierr);
164 
165   /* LU factorization - perms and factinfo are ignored by LAPACK */
166   i    = n-1;
167   ierr = MatZeroRows(mat,1,&i,-1.0,NULL,NULL);CHKERRQ(ierr);
168   ierr = MatMult(mat,x,b);CHKERRQ(ierr);
169 
170   /* in-place LU */
171   if (inplace) {
172     Mat RHS2;
173 
174     ierr = MatDuplicate(mat,MAT_COPY_VALUES,&F);CHKERRQ(ierr);
175     ierr = MatLUFactor(F,perm,perm,0);CHKERRQ(ierr);
176     ierr = MatSolve(F,b,y);CHKERRQ(ierr);
177     ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
178     ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
179     if (norm > tol) {
180       ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place LU %g\n",(double)norm);CHKERRQ(ierr);
181     }
182     ierr = MatMatSolve(F,RHS,SOLU);CHKERRQ(ierr);
183     ierr = MatMatMult(mat,SOLU,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RHS2);CHKERRQ(ierr);
184     ierr = MatAXPY(RHS,-1.0,RHS2,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
185     ierr = MatNorm(RHS,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
186     if (norm > tol) {
187       ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for in-place LU (MatMatSolve) %g\n",(double)norm);CHKERRQ(ierr);
188     }
189     ierr = MatDestroy(&F);CHKERRQ(ierr);
190     ierr = MatDestroy(&RHS2);CHKERRQ(ierr);
191   }
192 
193   /* out-of-place LU */
194   ierr = MatGetFactor(mat,solver,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
195   ierr = MatLUFactorSymbolic(F,mat,perm,perm,0);CHKERRQ(ierr);
196   ierr = MatLUFactorNumeric(F,mat,0);CHKERRQ(ierr);
197   ierr = MatSolve(F,b,y);CHKERRQ(ierr);
198   ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
199   ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
200   if (norm > tol) {
201     ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place LU %g\n",(double)norm);CHKERRQ(ierr);
202   }
203 
204   /* free space */
205   ierr = ISDestroy(&perm);CHKERRQ(ierr);
206   ierr = MatDestroy(&F);CHKERRQ(ierr);
207   ierr = MatDestroy(&mat);CHKERRQ(ierr);
208   ierr = MatDestroy(&RHS);CHKERRQ(ierr);
209   ierr = MatDestroy(&SOLU);CHKERRQ(ierr);
210   ierr = VecDestroy(&x);CHKERRQ(ierr);
211   ierr = VecDestroy(&b);CHKERRQ(ierr);
212   ierr = VecDestroy(&y);CHKERRQ(ierr);
213   ierr = VecDestroy(&ytmp);CHKERRQ(ierr);
214 
215   /* setup rectanglar */
216   ierr = createMatsAndVecs(m, n, nrhs, full, &mat, &RHS, &SOLU, &x, &y, &b);CHKERRQ(ierr);
217   ierr = VecDuplicate(y,&ytmp);CHKERRQ(ierr);
218 
219   /* QR factorization - perms and factinfo are ignored by LAPACK */
220   ierr = MatMult(mat,x,b);CHKERRQ(ierr);
221 
222   /* in-place QR */
223   if (inplace) {
224     Mat SOLU2;
225 
226     ierr = MatDuplicate(mat,MAT_COPY_VALUES,&F);CHKERRQ(ierr);
227     ierr = MatQRFactor(F,NULL,0);CHKERRQ(ierr);
228     ierr = MatSolve(F,b,y);CHKERRQ(ierr);
229     ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
230     ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
231     if (norm > tol) {
232       ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place QR %g\n",(double)norm);CHKERRQ(ierr);
233     }
234     ierr = MatMatMult(mat,SOLU,MAT_REUSE_MATRIX,PETSC_DEFAULT,&RHS);CHKERRQ(ierr);
235     ierr = MatDuplicate(SOLU, MAT_DO_NOT_COPY_VALUES, &SOLU2);CHKERRQ(ierr);
236     ierr = MatMatSolve(F,RHS,SOLU2);CHKERRQ(ierr);
237     ierr = MatAXPY(SOLU2,-1.0,SOLU,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
238     ierr = MatNorm(SOLU2,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
239     if (norm > tol) {
240       ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of error for in-place QR (MatMatSolve) %g\n",(double)norm);CHKERRQ(ierr);
241     }
242     ierr = MatDestroy(&F);CHKERRQ(ierr);
243     ierr = MatDestroy(&SOLU2);CHKERRQ(ierr);
244   }
245 
246   /* out-of-place QR */
247   ierr = MatGetFactor(mat,solver,MAT_FACTOR_QR,&F);CHKERRQ(ierr);
248   ierr = MatQRFactorSymbolic(F,mat,NULL,NULL);CHKERRQ(ierr);
249   ierr = MatQRFactorNumeric(F,mat,NULL);CHKERRQ(ierr);
250   ierr = MatSolve(F,b,y);CHKERRQ(ierr);
251   ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
252   ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
253   if (norm > tol) {
254     ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place QR %g\n",(double)norm);CHKERRQ(ierr);
255   }
256 
257   if (m == n) {
258     /* out-of-place MatSolveTranspose */
259     ierr = MatMultTranspose(mat,x,b);CHKERRQ(ierr);
260     ierr = MatSolveTranspose(F,b,y);CHKERRQ(ierr);
261     ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr);
262     ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
263     if (norm > tol) {
264       ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place QR %g\n",(double)norm);CHKERRQ(ierr);
265     }
266   }
267 
268   /* free space */
269   ierr = MatDestroy(&F);CHKERRQ(ierr);
270   ierr = MatDestroy(&mat);CHKERRQ(ierr);
271   ierr = MatDestroy(&RHS);CHKERRQ(ierr);
272   ierr = MatDestroy(&SOLU);CHKERRQ(ierr);
273   ierr = VecDestroy(&x);CHKERRQ(ierr);
274   ierr = VecDestroy(&b);CHKERRQ(ierr);
275   ierr = VecDestroy(&y);CHKERRQ(ierr);
276   ierr = VecDestroy(&ytmp);CHKERRQ(ierr);
277   ierr = PetscFinalize();
278   return ierr;
279 }
280 
281 
282 
283 /*TEST
284 
285    test:
286 
287    test:
288      requires: cuda
289      suffix: seqdensecuda
290      args: -mat_type seqdensecuda -rhs_mat_type seqdensecuda -ldl 0 -solver_type {{petsc cuda}}
291      output_file: output/ex1_1.out
292 
293    test:
294      requires: cuda
295      suffix: seqdensecuda_seqaijcusparse
296      args: -mat_type seqaijcusparse -rhs_mat_type seqdensecuda
297      output_file: output/ex1_2.out
298 
299    test:
300      requires: cuda viennacl
301      suffix: seqdensecuda_seqaijviennacl
302      args: -mat_type seqaijviennacl -rhs_mat_type seqdensecuda
303      output_file: output/ex1_2.out
304 
305    test:
306      suffix: 4
307      args: -m 10 -n 10
308      output_file: output/ex1_1.out
309 
310 TEST*/
311