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_DEFAULT, &C));
127:     } else {
128:       Mat Qc;
129:       PetscCall(MatHermitianTranspose(Q, MAT_INITIAL_MATRIX, &Qc));
130:       PetscCall(MatMatMult(Qc, Q, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &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*/