1 static char help[] = "Tests MatSolve(), MatSolveTranspose() and MatMatSolve() with SEQDENSE\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc,char **args)
6 {
7 Mat A,RHS,C,F,X;
8 Vec u,x,b;
9 PetscErrorCode ierr;
10 PetscMPIInt size;
11 PetscInt m,n,nsolve,nrhs;
12 PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON;
13 PetscRandom rand;
14 PetscBool data_provided,herm,symm,hpd;
15 MatFactorType ftyp;
16 PetscViewer fd;
17 char file[PETSC_MAX_PATH_LEN];
18
19 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
21 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor test");
22 /* Determine which type of solver we want to test for */
23 herm = PETSC_FALSE;
24 symm = PETSC_FALSE;
25 hpd = PETSC_FALSE;
26 ierr = PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL);CHKERRQ(ierr);
27 ierr = PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL);CHKERRQ(ierr);
28 ierr = PetscOptionsGetBool(NULL,NULL,"-hpd_solve",&hpd,NULL);CHKERRQ(ierr);
29
30 /* Determine file from which we read the matrix A */
31 ftyp = MAT_FACTOR_LU;
32 ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&data_provided);CHKERRQ(ierr);
33 if (!data_provided) { /* get matrices from PETSc distribution */
34 ierr = PetscStrcpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/");CHKERRQ(ierr);
35 if (hpd) {
36 #if defined(PETSC_USE_COMPLEX)
37 ierr = PetscStrcat(file,"hpd-complex-");CHKERRQ(ierr);
38 #else
39 ierr = PetscStrcat(file,"spd-real-");CHKERRQ(ierr);
40 #endif
41 ftyp = MAT_FACTOR_CHOLESKY;
42 } else {
43 #if defined(PETSC_USE_COMPLEX)
44 ierr = PetscStrcat(file,"nh-complex-");CHKERRQ(ierr);
45 #else
46 ierr = PetscStrcat(file,"ns-real-");CHKERRQ(ierr);
47 #endif
48 }
49 #if defined(PETSC_USE_64BIT_INDICES)
50 ierr = PetscStrcat(file,"int64-");CHKERRQ(ierr);
51 #else
52 ierr = PetscStrcat(file,"int32-");CHKERRQ(ierr);
53 #endif
54 #if defined(PETSC_USE_REAL_SINGLE)
55 ierr = PetscStrcat(file,"float32");CHKERRQ(ierr);
56 #else
57 ierr = PetscStrcat(file,"float64");CHKERRQ(ierr);
58 #endif
59 }
60
61 /* Load matrix A */
62 #if defined(PETSC_USE_REAL___FLOAT128)
63 ierr = PetscOptionsInsertString(NULL,"-binary_read_double");CHKERRQ(ierr);
64 #endif
65 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
66 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
67 ierr = MatLoad(A,fd);CHKERRQ(ierr);
68 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
69 ierr = MatConvert(A,MATSEQDENSE,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr);
70 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
71 if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
72
73 /* Create dense matrix C and X; C holds true solution with identical colums */
74 nrhs = 2;
75 ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr);
76 ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
77 ierr = MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);CHKERRQ(ierr);
78 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr);
79 ierr = MatSetFromOptions(C);CHKERRQ(ierr);
80 ierr = MatSetUp(C);CHKERRQ(ierr);
81
82 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr);
83 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr);
84 ierr = MatSetRandom(C,rand);CHKERRQ(ierr);
85 ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);CHKERRQ(ierr);
86 ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&RHS);CHKERRQ(ierr);
87
88 /* Create vectors */
89 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
90 ierr = VecSetSizes(x,n,PETSC_DECIDE);CHKERRQ(ierr);
91 ierr = VecSetFromOptions(x);CHKERRQ(ierr);
92 ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
93 ierr = VecDuplicate(x,&u);CHKERRQ(ierr); /* save the true solution */
94
95 /* make a symmetric matrix */
96 if (symm) {
97 Mat AT;
98
99 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&AT);CHKERRQ(ierr);
100 ierr = MatAXPY(A,1.0,AT,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
101 ierr = MatDestroy(&AT);CHKERRQ(ierr);
102 ftyp = MAT_FACTOR_CHOLESKY;
103 }
104 /* make an hermitian matrix */
105 if (herm) {
106 Mat AH;
107
108 ierr = MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&AH);CHKERRQ(ierr);
109 ierr = MatAXPY(A,1.0,AH,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
110 ierr = MatDestroy(&AH);CHKERRQ(ierr);
111 ftyp = MAT_FACTOR_CHOLESKY;
112 }
113 ierr = PetscObjectSetName((PetscObject)A,"A");CHKERRQ(ierr);
114 ierr = MatViewFromOptions(A,NULL,"-amat_view");CHKERRQ(ierr);
115
116 ierr = MatDuplicate(A,MAT_COPY_VALUES,&F);CHKERRQ(ierr);
117 ierr = MatSetOption(F,MAT_SYMMETRIC,symm);CHKERRQ(ierr);
118 /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */
119 ierr = MatSetOption(F,MAT_HERMITIAN,(PetscBool)(hpd || herm));CHKERRQ(ierr);
120 ierr = MatSetOption(F,MAT_SPD,hpd);CHKERRQ(ierr);
121 if (ftyp == MAT_FACTOR_LU) {
122 ierr = MatLUFactor(F,NULL,NULL,NULL);CHKERRQ(ierr);
123 } else {
124 ierr = MatCholeskyFactor(F,NULL,NULL);CHKERRQ(ierr);
125 }
126
127 for (nsolve = 0; nsolve < 2; nsolve++) {
128 ierr = VecSetRandom(x,rand);CHKERRQ(ierr);
129 ierr = VecCopy(x,u);CHKERRQ(ierr);
130 if (nsolve) {
131 ierr = MatMult(A,x,b);CHKERRQ(ierr);
132 ierr = MatSolve(F,b,x);CHKERRQ(ierr);
133 } else {
134 ierr = MatMultTranspose(A,x,b);CHKERRQ(ierr);
135 ierr = MatSolveTranspose(F,b,x);CHKERRQ(ierr);
136 }
137 /* Check the error */
138 ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr); /* u <- (-1.0)x + u */
139 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
140 if (norm > tol) {
141 PetscReal resi;
142 if (nsolve) {
143 ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */
144 } else {
145 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); /* u = A*x */
146 }
147 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); /* u <- (-1.0)b + u */
148 ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr);
149 if (nsolve) {
150 ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolve error: Norm of error %g, residual %f\n",norm,resi);CHKERRQ(ierr);
151 } else {
152 ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolveTranspose error: Norm of error %g, residual %f\n",norm,resi);CHKERRQ(ierr);
153 }
154 }
155 }
156 ierr = MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr);
157 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr);
158
159 /* Check the error */
160 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
161 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
162 if (norm > tol) {
163 ierr = PetscPrintf(PETSC_COMM_SELF,"MatMatSolve: Norm of error %g\n",norm);CHKERRQ(ierr);
164 }
165
166 /* Free data structures */
167 ierr = MatDestroy(&A);CHKERRQ(ierr);
168 ierr = MatDestroy(&C);CHKERRQ(ierr);
169 ierr = MatDestroy(&F);CHKERRQ(ierr);
170 ierr = MatDestroy(&X);CHKERRQ(ierr);
171 ierr = MatDestroy(&RHS);CHKERRQ(ierr);
172 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr);
173 ierr = VecDestroy(&x);CHKERRQ(ierr);
174 ierr = VecDestroy(&b);CHKERRQ(ierr);
175 ierr = VecDestroy(&u);CHKERRQ(ierr);
176 ierr = PetscFinalize();
177 return ierr;
178 }
179
180
181 /*TEST
182
183 testset:
184 output_file: output/ex215.out
185 test:
186 suffix: ns
187 test:
188 suffix: sym
189 args: -symmetric_solve
190 test:
191 suffix: herm
192 args: -hermitian_solve
193 test:
194 suffix: hpd
195 args: -hpd_solve
196
197 TEST*/
198