Actual source code: ex54.c
1: /*
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: */
7: #include <petscksp.h>
9: PetscErrorCode fill(Mat m, Vec v)
10: {
11: PetscInt idxn[3] = {0, 1, 2};
12: PetscInt localRows = 0;
13: PetscMPIInt rank, size;
15: PetscFunctionBegin;
16: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
17: PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
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));
24: PetscCall(MatSetFromOptions(m));
25: PetscCall(VecSetFromOptions(v));
26: PetscCall(MatSetUp(m));
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};
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));
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};
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};
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: }
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;
81: PetscFunctionBeginUser;
82: PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
84: PetscCall(VecCreate(PETSC_COMM_WORLD, &v));
85: PetscCall(MatCreate(PETSC_COMM_WORLD, &Q));
86: PetscCall(MatSetType(Q, MATDENSE));
87: PetscCall(fill(Q, v));
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));
156: PetscCall(PetscFinalize());
157: return 0;
158: }
160: /*TEST
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
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
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
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
180: test:
181: requires: suitesparse
182: suffix: 5
183: nsize: 1
184: args: -exact
186: TEST*/