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*/