1 /*
2
3 Tests MPIDENSE matrix operations MatMultTranspose() with processes with no rows or columns.
4 As the matrix is rectangular, least square solution is computed, so KSPLSQR is also tested here.
5 */
6
7 #include <petscksp.h>
8
fill(Mat m,Vec v)9 PetscErrorCode fill(Mat m, Vec v)
10 {
11 PetscInt idxn[3] = {0, 1, 2};
12 PetscInt localRows = 0;
13 PetscMPIInt rank, size;
14
15 PetscFunctionBegin;
16 PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
17 PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
18
19 if (rank == 1 || rank == 2) localRows = 4;
20 if (size == 1) localRows = 8;
21 PetscCall(MatSetSizes(m, localRows, PETSC_DECIDE, PETSC_DECIDE, 3));
22 PetscCall(VecSetSizes(v, localRows, PETSC_DECIDE));
23
24 PetscCall(MatSetFromOptions(m));
25 PetscCall(VecSetFromOptions(v));
26 PetscCall(MatSetUp(m));
27
28 if (size == 1) {
29 PetscInt idxm1[4] = {0, 1, 2, 3};
30 PetscScalar values1[12] = {1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1};
31 PetscInt idxm2[4] = {4, 5, 6, 7};
32 PetscScalar values2[12] = {1, 2, 0, 1, 2, 1, 1, 3, 0, 1, 3, 1};
33
34 PetscCall(MatSetValues(m, 4, idxm1, 3, idxn, values1, INSERT_VALUES));
35 PetscCall(VecSetValue(v, 0, 1.1, INSERT_VALUES));
36 PetscCall(VecSetValue(v, 1, 2.5, INSERT_VALUES));
37 PetscCall(VecSetValue(v, 2, 3, INSERT_VALUES));
38 PetscCall(VecSetValue(v, 3, 4, INSERT_VALUES));
39
40 PetscCall(MatSetValues(m, 4, idxm2, 3, idxn, values2, INSERT_VALUES));
41 PetscCall(VecSetValue(v, 4, 5, INSERT_VALUES));
42 PetscCall(VecSetValue(v, 5, 6, INSERT_VALUES));
43 PetscCall(VecSetValue(v, 6, 7, INSERT_VALUES));
44 PetscCall(VecSetValue(v, 7, 8, INSERT_VALUES));
45 } else if (rank == 1) {
46 PetscInt idxm[4] = {0, 1, 2, 3};
47 PetscScalar values[12] = {1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1};
48
49 PetscCall(MatSetValues(m, 4, idxm, 3, idxn, values, INSERT_VALUES));
50 PetscCall(VecSetValue(v, 0, 1.1, INSERT_VALUES));
51 PetscCall(VecSetValue(v, 1, 2.5, INSERT_VALUES));
52 PetscCall(VecSetValue(v, 2, 3, INSERT_VALUES));
53 PetscCall(VecSetValue(v, 3, 4, INSERT_VALUES));
54 } else if (rank == 2) {
55 PetscInt idxm[4] = {4, 5, 6, 7};
56 PetscScalar values[12] = {1, 2, 0, 1, 2, 1, 1, 3, 0, 1, 3, 1};
57
58 PetscCall(MatSetValues(m, 4, idxm, 3, idxn, values, INSERT_VALUES));
59 PetscCall(VecSetValue(v, 4, 5, INSERT_VALUES));
60 PetscCall(VecSetValue(v, 5, 6, INSERT_VALUES));
61 PetscCall(VecSetValue(v, 6, 7, INSERT_VALUES));
62 PetscCall(VecSetValue(v, 7, 8, INSERT_VALUES));
63 }
64 PetscCall(MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY));
65 PetscCall(MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY));
66 PetscCall(VecAssemblyBegin(v));
67 PetscCall(VecAssemblyEnd(v));
68 PetscFunctionReturn(PETSC_SUCCESS);
69 }
70
main(int argc,char ** argv)71 int main(int argc, char **argv)
72 {
73 Mat Q, C, V, A, B;
74 Vec v, a, b, se;
75 KSP QRsolver;
76 PC pc;
77 PetscReal norm;
78 PetscInt m, n;
79 PetscBool exact = PETSC_FALSE;
80
81 PetscFunctionBeginUser;
82 PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
83
84 PetscCall(VecCreate(PETSC_COMM_WORLD, &v));
85 PetscCall(MatCreate(PETSC_COMM_WORLD, &Q));
86 PetscCall(MatSetType(Q, MATDENSE));
87 PetscCall(fill(Q, v));
88
89 PetscCall(MatCreateVecs(Q, &a, NULL));
90 PetscCall(MatCreateNormalHermitian(Q, &C));
91 PetscCall(KSPCreate(PETSC_COMM_WORLD, &QRsolver));
92 PetscCall(KSPGetPC(QRsolver, &pc));
93 PetscCall(PCSetType(pc, PCNONE));
94 PetscCall(KSPSetType(QRsolver, KSPLSQR));
95 PetscCall(KSPSetFromOptions(QRsolver));
96 PetscCall(KSPSetOperators(QRsolver, Q, Q));
97 PetscCall(MatViewFromOptions(Q, NULL, "-sys_view"));
98 PetscCall(VecViewFromOptions(a, NULL, "-rhs_view"));
99 PetscCall(KSPSolve(QRsolver, v, a));
100 PetscCall(KSPLSQRGetStandardErrorVec(QRsolver, &se));
101 if (se) PetscCall(VecViewFromOptions(se, NULL, "-se_view"));
102 PetscCall(PetscOptionsGetBool(NULL, NULL, "-exact", &exact, NULL));
103 if (exact) {
104 PetscCall(KSPDestroy(&QRsolver));
105 PetscCall(MatDestroy(&C));
106 PetscCall(MatConvert(Q, MATAIJ, MAT_INPLACE_MATRIX, &Q));
107 PetscCall(MatCreateNormalHermitian(Q, &C));
108 PetscCall(KSPCreate(PETSC_COMM_WORLD, &QRsolver));
109 PetscCall(KSPGetPC(QRsolver, &pc));
110 PetscCall(PCSetType(pc, PCQR));
111 PetscCall(KSPSetType(QRsolver, KSPLSQR));
112 PetscCall(KSPSetFromOptions(QRsolver));
113 PetscCall(KSPSetOperators(QRsolver, Q, C));
114 PetscCall(VecZeroEntries(a));
115 PetscCall(KSPSolve(QRsolver, v, a));
116 PetscCall(MatGetLocalSize(Q, &m, &n));
117 PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, 5, NULL, &V));
118 PetscCall(MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DECIDE, 5, NULL, &A));
119 PetscCall(MatDuplicate(A, MAT_SHARE_NONZERO_PATTERN, &B));
120 PetscCall(MatSetRandom(V, NULL));
121 PetscCall(KSPMatSolve(QRsolver, V, A));
122 PetscCall(KSPView(QRsolver, PETSC_VIEWER_STDOUT_WORLD));
123 PetscCall(PCSetType(pc, PCCHOLESKY));
124 PetscCall(MatDestroy(&C));
125 if (!PetscDefined(USE_COMPLEX)) {
126 PetscCall(MatTransposeMatMult(Q, Q, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
127 } else {
128 Mat Qc;
129 PetscCall(MatHermitianTranspose(Q, MAT_INITIAL_MATRIX, &Qc));
130 PetscCall(MatMatMult(Qc, Q, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
131 PetscCall(MatDestroy(&Qc));
132 }
133 PetscCall(KSPSetOperators(QRsolver, Q, C));
134 PetscCall(KSPSetFromOptions(QRsolver));
135 PetscCall(VecDuplicate(a, &b));
136 PetscCall(KSPSolve(QRsolver, v, b));
137 PetscCall(KSPMatSolve(QRsolver, V, B));
138 PetscCall(KSPView(QRsolver, PETSC_VIEWER_STDOUT_WORLD));
139 PetscCall(VecAXPY(a, -1.0, b));
140 PetscCall(VecNorm(a, NORM_2, &norm));
141 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||a-b|| > PETSC_SMALL (%g)", (double)norm);
142 PetscCall(MatAXPY(A, -1.0, B, SAME_NONZERO_PATTERN));
143 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm));
144 PetscCheck(norm <= PETSC_SMALL, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "||A-B|| > PETSC_SMALL (%g)", (double)norm);
145 PetscCall(VecDestroy(&b));
146 PetscCall(MatDestroy(&V));
147 PetscCall(MatDestroy(&A));
148 PetscCall(MatDestroy(&B));
149 }
150 PetscCall(KSPDestroy(&QRsolver));
151 PetscCall(VecDestroy(&a));
152 PetscCall(VecDestroy(&v));
153 PetscCall(MatDestroy(&C));
154 PetscCall(MatDestroy(&Q));
155
156 PetscCall(PetscFinalize());
157 return 0;
158 }
159
160 /*TEST
161
162 test:
163 args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_compute_standard_error -ksp_lsqr_monitor
164
165 test:
166 suffix: 2
167 nsize: 4
168 args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_compute_standard_error -ksp_lsqr_monitor
169
170 test:
171 suffix: 3
172 nsize: 2
173 args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_monitor -ksp_convergence_test lsqr -ksp_lsqr_compute_standard_error -se_view -ksp_lsqr_exact_mat_norm 0
174
175 test:
176 suffix: 4
177 nsize: 2
178 args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_monitor -ksp_convergence_test lsqr -ksp_lsqr_compute_standard_error -se_view -ksp_lsqr_exact_mat_norm 1
179
180 test:
181 requires: suitesparse
182 suffix: 5
183 nsize: 1
184 args: -exact
185
186 TEST*/
187