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: MPI_Comm_rank(MPI_COMM_WORLD, &rank);
16: MPI_Comm_size(MPI_COMM_WORLD, &size);
18: if (rank == 1 || rank == 2) localRows = 4;
19: if (size == 1) localRows = 8;
20: MatSetSizes(m, localRows, PETSC_DECIDE, PETSC_DECIDE, 3);
21: VecSetSizes(v, localRows, PETSC_DECIDE);
23: MatSetFromOptions(m);
24: VecSetFromOptions(v);
25: MatSetUp(m);
27: if (size == 1) {
28: PetscInt idxm1[4] = {0, 1, 2, 3};
29: PetscScalar values1[12] = {1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1};
30: PetscInt idxm2[4] = {4, 5, 6, 7};
31: PetscScalar values2[12] = {1, 2, 0, 1, 2, 1, 1, 3, 0, 1, 3, 1};
33: MatSetValues(m, 4, idxm1, 3, idxn, values1, INSERT_VALUES);
34: VecSetValue(v, 0, 1.1, INSERT_VALUES); VecSetValue(v, 1, 2.5, INSERT_VALUES);
35: VecSetValue(v, 2, 3, INSERT_VALUES); VecSetValue(v, 3, 4, INSERT_VALUES);
37: MatSetValues(m, 4, idxm2, 3, idxn, values2, INSERT_VALUES);
38: VecSetValue(v, 4, 5, INSERT_VALUES); VecSetValue(v, 5, 6, INSERT_VALUES);
39: VecSetValue(v, 6, 7, INSERT_VALUES); VecSetValue(v, 7, 8, INSERT_VALUES);
40: } else if (rank == 1) {
41: PetscInt idxm[4] = {0, 1, 2, 3};
42: PetscScalar values[12] = {1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1};
44: MatSetValues(m, 4, idxm, 3, idxn, values, INSERT_VALUES);
45: VecSetValue(v, 0, 1.1, INSERT_VALUES); VecSetValue(v, 1, 2.5, INSERT_VALUES);
46: VecSetValue(v, 2, 3, INSERT_VALUES); VecSetValue(v, 3, 4, INSERT_VALUES);
47: } else if (rank == 2) {
48: PetscInt idxm[4] = {4, 5, 6, 7};
49: PetscScalar values[12] = {1, 2, 0, 1, 2, 1, 1, 3, 0, 1, 3, 1};
51: MatSetValues(m, 4, idxm, 3, idxn, values, INSERT_VALUES);
52: VecSetValue(v, 4, 5, INSERT_VALUES); VecSetValue(v, 5, 6, INSERT_VALUES);
53: VecSetValue(v, 6, 7, INSERT_VALUES); VecSetValue(v, 7, 8, INSERT_VALUES);
54: }
55: MatAssemblyBegin(m, MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(m, MAT_FINAL_ASSEMBLY);
57: VecAssemblyBegin(v);
58: VecAssemblyEnd(v);
59: return 0;
60: }
62: int main(int argc, char **argv)
63: {
64: Mat Q, C, V, A, B;
65: Vec v, a, b, se;
66: KSP QRsolver;
67: PC pc;
68: PetscReal norm;
69: PetscInt m, n;
70: PetscBool exact = PETSC_FALSE;
73: PetscInitialize(&argc, &argv, NULL, NULL);
75: VecCreate(PETSC_COMM_WORLD, &v);
76: MatCreate(PETSC_COMM_WORLD, &Q);
77: MatSetType(Q, MATDENSE);
78: fill(Q, v);
80: MatCreateVecs(Q, &a, NULL);
81: MatCreateNormalHermitian(Q, &C);
82: KSPCreate(PETSC_COMM_WORLD, &QRsolver);
83: KSPGetPC(QRsolver, &pc);
84: PCSetType(pc, PCNONE);
85: KSPSetType(QRsolver, KSPLSQR);
86: KSPSetFromOptions(QRsolver);
87: KSPSetOperators(QRsolver, Q, Q);
88: MatViewFromOptions(Q, NULL, "-sys_view");
89: VecViewFromOptions(a, NULL, "-rhs_view");
90: KSPSolve(QRsolver, v, a);
91: KSPLSQRGetStandardErrorVec(QRsolver, &se);
92: if (se) VecViewFromOptions(se, NULL, "-se_view");
93: PetscOptionsGetBool(NULL, NULL, "-exact", &exact, NULL);
94: if (exact) {
95: KSPDestroy(&QRsolver);
96: MatDestroy(&C);
97: MatConvert(Q, MATAIJ, MAT_INPLACE_MATRIX, &Q);
98: MatCreateNormalHermitian(Q, &C);
99: KSPCreate(PETSC_COMM_WORLD, &QRsolver);
100: KSPGetPC(QRsolver, &pc);
101: PCSetType(pc, PCQR);
102: KSPSetType(QRsolver, KSPLSQR);
103: KSPSetFromOptions(QRsolver);
104: KSPSetOperators(QRsolver, Q, C);
105: VecZeroEntries(a);
106: KSPSolve(QRsolver, v, a);
107: MatGetLocalSize(Q, &m, &n);
108: MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, 5, NULL, &V);
109: MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DECIDE, 5, NULL, &A);
110: MatDuplicate(A, MAT_SHARE_NONZERO_PATTERN, &B);
111: MatSetRandom(V, NULL);
112: KSPMatSolve(QRsolver, V, A);
113: KSPView(QRsolver, PETSC_VIEWER_STDOUT_WORLD);
114: PCSetType(pc, PCCHOLESKY);
115: MatDestroy(&C);
116: if (!PetscDefined(USE_COMPLEX)) {
117: MatTransposeMatMult(Q, Q, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
118: } else {
119: Mat Qc;
120: MatHermitianTranspose(Q, MAT_INITIAL_MATRIX, &Qc);
121: MatMatMult(Qc, Q, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
122: MatDestroy(&Qc);
123: }
124: KSPSetOperators(QRsolver, Q, C);
125: KSPSetFromOptions(QRsolver);
126: VecDuplicate(a, &b);
127: KSPSolve(QRsolver, v, b);
128: KSPMatSolve(QRsolver, V, B);
129: KSPView(QRsolver, PETSC_VIEWER_STDOUT_WORLD);
130: VecAXPY(a, -1.0, b);
131: VecNorm(a, NORM_2, &norm);
133: MatAXPY(A, -1.0, B, SAME_NONZERO_PATTERN);
134: MatNorm(A, NORM_FROBENIUS, &norm);
136: VecDestroy(&b);
137: MatDestroy(&V);
138: MatDestroy(&A);
139: MatDestroy(&B);
140: }
141: KSPDestroy(&QRsolver);
142: VecDestroy(&a);
143: VecDestroy(&v);
144: MatDestroy(&C);
145: MatDestroy(&Q);
147: PetscFinalize();
148: return 0;
149: }
151: /*TEST
153: test:
154: args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_compute_standard_error -ksp_lsqr_monitor
156: test:
157: suffix: 2
158: nsize: 4
159: args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_compute_standard_error -ksp_lsqr_monitor
161: test:
162: suffix: 3
163: nsize: 2
164: 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
166: test:
167: suffix: 4
168: nsize: 2
169: 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
171: test:
172: requires: suitesparse
173: suffix: 5
174: nsize: 1
175: args: -exact
177: TEST*/