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