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;
72: PetscInitialize(&argc, &argv, NULL, NULL);
74: VecCreate(PETSC_COMM_WORLD, &v);
75: MatCreate(PETSC_COMM_WORLD, &Q);
76: MatSetType(Q, MATDENSE);
77: fill(Q, v);
79: MatCreateVecs(Q, &a, NULL);
80: MatCreateNormalHermitian(Q, &C);
81: KSPCreate(PETSC_COMM_WORLD, &QRsolver);
82: KSPGetPC(QRsolver, &pc);
83: PCSetType(pc, PCNONE);
84: KSPSetType(QRsolver, KSPLSQR);
85: KSPSetFromOptions(QRsolver);
86: KSPSetOperators(QRsolver, Q, Q);
87: MatViewFromOptions(Q, NULL, "-sys_view");
88: VecViewFromOptions(a, NULL, "-rhs_view");
89: KSPSolve(QRsolver, v, a);
90: KSPLSQRGetStandardErrorVec(QRsolver, &se);
91: if (se) {
92: VecViewFromOptions(se, NULL, "-se_view");
93: }
94: PetscOptionsGetBool(NULL, NULL, "-exact", &exact, NULL);
95: if (exact) {
96: KSPDestroy(&QRsolver);
97: MatDestroy(&C);
98: MatConvert(Q, MATAIJ, MAT_INPLACE_MATRIX, &Q);
99: MatCreateNormalHermitian(Q, &C);
100: KSPCreate(PETSC_COMM_WORLD, &QRsolver);
101: KSPGetPC(QRsolver, &pc);
102: PCSetType(pc, PCQR);
103: KSPSetType(QRsolver, KSPLSQR);
104: KSPSetFromOptions(QRsolver);
105: KSPSetOperators(QRsolver, Q, C);
106: VecZeroEntries(a);
107: KSPSolve(QRsolver, v, a);
108: MatGetLocalSize(Q, &m, &n);
109: MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, 5, NULL, &V);
110: MatCreateDense(PETSC_COMM_WORLD, n, PETSC_DECIDE, PETSC_DECIDE, 5, NULL, &A);
111: MatDuplicate(A, MAT_SHARE_NONZERO_PATTERN, &B);
112: MatSetRandom(V, NULL);
113: KSPMatSolve(QRsolver, V, A);
114: KSPView(QRsolver, PETSC_VIEWER_STDOUT_WORLD);
115: PCSetType(pc, PCCHOLESKY);
116: MatDestroy(&C);
117: if (!PetscDefined(USE_COMPLEX)) {
118: MatTransposeMatMult(Q, Q, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
119: } else {
120: Mat Qc;
121: MatHermitianTranspose(Q, MAT_INITIAL_MATRIX, &Qc);
122: MatMatMult(Qc, Q, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
123: MatDestroy(&Qc);
124: }
125: KSPSetOperators(QRsolver, Q, C);
126: KSPSetFromOptions(QRsolver);
127: VecDuplicate(a, &b);
128: KSPSolve(QRsolver, v, b);
129: KSPMatSolve(QRsolver, V, B);
130: KSPView(QRsolver, PETSC_VIEWER_STDOUT_WORLD);
131: VecAXPY(a, -1.0, b);
132: VecNorm(a, NORM_2, &norm);
134: MatAXPY(A, -1.0, B, SAME_NONZERO_PATTERN);
135: MatNorm(A, NORM_FROBENIUS, &norm);
137: VecDestroy(&b);
138: MatDestroy(&V);
139: MatDestroy(&A);
140: MatDestroy(&B);
141: }
142: KSPDestroy(&QRsolver);
143: VecDestroy(&a);
144: VecDestroy(&v);
145: MatDestroy(&C);
146: MatDestroy(&Q);
148: PetscFinalize();
149: return 0;
150: }
152: /*TEST
154: test:
155: args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_compute_standard_error -ksp_lsqr_monitor
157: test:
158: suffix: 2
159: nsize: 4
160: args: -ksp_monitor_true_residual -ksp_max_it 10 -sys_view -ksp_converged_reason -ksp_view -ksp_lsqr_compute_standard_error -ksp_lsqr_monitor
162: test:
163: suffix: 3
164: nsize: 2
165: 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
167: test:
168: suffix: 4
169: nsize: 2
170: 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
172: test:
173: requires: suitesparse
174: suffix: 5
175: nsize: 1
176: args: -exact
178: TEST*/