Actual source code: lmvm_copy_test.c

  1: const char help[] = "Test that MatCopy() does not affect the copied LMVM matrix";

  3: #include <petsc.h>

  5: static PetscErrorCode positiveVectorUpdate(PetscRandom rand, Vec x, Vec f)
  6: {
  7:   Vec         _x, _f;
  8:   PetscScalar dot;

 10:   PetscFunctionBegin;
 11:   PetscCall(VecDuplicate(x, &_x));
 12:   PetscCall(VecDuplicate(f, &_f));
 13:   PetscCall(VecSetRandom(_x, rand));
 14:   PetscCall(VecSetRandom(_f, rand));
 15:   PetscCall(VecDot(_x, _f, &dot));
 16:   PetscCall(VecAXPY(x, PetscAbsScalar(dot) / dot, _x));
 17:   PetscCall(VecAXPY(f, 1.0, _f));
 18:   PetscCall(VecDestroy(&_f));
 19:   PetscCall(VecDestroy(&_x));
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

 23: static PetscErrorCode VecEqualToTolerance(Vec a, Vec b, NormType norm_type, PetscReal tol, PetscBool *flg)
 24: {
 25:   Vec       diff;
 26:   PetscReal diff_norm;

 28:   PetscFunctionBegin;
 29:   PetscCall(VecDuplicate(a, &diff));
 30:   PetscCall(VecCopy(a, diff));
 31:   PetscCall(VecAXPY(diff, -1.0, b));
 32:   PetscCall(VecNorm(diff, norm_type, &diff_norm));
 33:   PetscCall(VecDestroy(&diff));
 34:   *flg = (diff_norm <= tol) ? PETSC_TRUE : PETSC_FALSE;
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

 38: // unlike MatTestEqual(), this test tests MatMult() and MatSolve()
 39: static PetscErrorCode testMatEqual(PetscRandom rand, Mat A, Mat B, PetscBool *flg)
 40: {
 41:   Vec x, y_A, y_B;

 43:   PetscFunctionBegin;
 44:   *flg = PETSC_TRUE;
 45:   PetscCall(MatCreateVecs(A, &x, &y_A));
 46:   PetscCall(MatCreateVecs(B, NULL, &y_B));
 47:   PetscCall(VecSetRandom(x, rand));
 48:   PetscCall(MatMult(A, x, y_A));
 49:   PetscCall(MatMult(B, x, y_B));
 50:   PetscCall(VecEqualToTolerance(y_A, y_B, NORM_2, PETSC_SMALL, flg));
 51:   if (*flg == PETSC_TRUE) {
 52:     PetscCall(MatSolve(A, x, y_A));
 53:     PetscCall(MatSolve(B, x, y_B));
 54:     PetscCall(VecEqualToTolerance(y_A, y_B, NORM_2, PETSC_SMALL, flg));
 55:     if (*flg == PETSC_FALSE) {
 56:       PetscReal norm;

 58:       PetscCall(VecAXPY(y_A, -1.0, y_B));
 59:       PetscCall(VecNorm(y_A, NORM_INFINITY, &norm));
 60:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)A), "MatSolve() norm error %g\n", (double)norm));
 61:     }
 62:   } else {
 63:     PetscReal norm;

 65:     PetscCall(VecAXPY(y_A, -1.0, y_B));
 66:     PetscCall(VecNorm(y_A, NORM_INFINITY, &norm));
 67:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)A), "MatMult() norm error %g\n", (double)norm));
 68:   }
 69:   PetscCall(VecDestroy(&y_B));
 70:   PetscCall(VecDestroy(&y_A));
 71:   PetscCall(VecDestroy(&x));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: static PetscErrorCode testUnchangedBegin(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z)
 76: {
 77:   Vec _x, _y, _z;

 79:   PetscFunctionBegin;
 80:   PetscCall(MatCreateVecs(A, &_x, &_y));
 81:   PetscCall(MatCreateVecs(A, NULL, &_z));
 82:   PetscCall(VecSetRandom(_x, rand));
 83:   PetscCall(MatMult(A, _x, _y));
 84:   PetscCall(MatSolve(A, _x, _z));
 85:   *x = _x;
 86:   *y = _y;
 87:   *z = _z;
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: static PetscErrorCode testUnchangedEnd(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z, PetscBool *unchanged)
 92: {
 93:   Vec _x, _y, _z, _y2, _z2;

 95:   PetscFunctionBegin;
 96:   *unchanged = PETSC_TRUE;
 97:   _x         = *x;
 98:   _y         = *y;
 99:   _z         = *z;
100:   *x         = NULL;
101:   *y         = NULL;
102:   *z         = NULL;
103:   PetscCall(MatCreateVecs(A, NULL, &_y2));
104:   PetscCall(MatCreateVecs(A, NULL, &_z2));
105:   PetscCall(MatMult(A, _x, _y2));
106:   PetscCall(MatSolve(A, _x, _z2));
107:   PetscCall(VecEqual(_y, _y2, unchanged));
108:   if (*unchanged == PETSC_TRUE) PetscCall(VecEqual(_z, _z2, unchanged));
109:   PetscCall(VecDestroy(&_z2));
110:   PetscCall(VecDestroy(&_y2));
111:   PetscCall(VecDestroy(&_z));
112:   PetscCall(VecDestroy(&_y));
113:   PetscCall(VecDestroy(&_x));
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static PetscErrorCode testMatLMVMCopy(PetscRandom rand)
118: {
119:   PetscInt    N     = 10;
120:   MPI_Comm    comm  = PetscObjectComm((PetscObject)rand);
121:   PetscInt    k_pre = 2; // number of updates before copy
122:   Mat         A, A_copy;
123:   Vec         u, x, f, x_copy, f_copy, v1, v2, v3;
124:   PetscBool   equal;
125:   PetscLayout layout;

127:   PetscFunctionBegin;
128:   PetscCall(VecCreateMPI(comm, PETSC_DECIDE, N, &u));
129:   PetscCall(VecSetFromOptions(u));
130:   PetscCall(VecSetUp(u));
131:   PetscCall(VecDuplicate(u, &x));
132:   PetscCall(VecDuplicate(u, &f));
133:   PetscCall(VecGetLayout(u, &layout));
134:   PetscCall(MatCreate(comm, &A));
135:   PetscCall(MatSetLayouts(A, layout, layout));
136:   PetscCall(MatSetType(A, MATLMVMBFGS));
137:   PetscCall(MatSetFromOptions(A));
138:   PetscCall(positiveVectorUpdate(rand, x, f));
139:   PetscCall(MatLMVMAllocate(A, x, f));
140:   PetscCall(MatSetUp(A));
141:   for (PetscInt k = 0; k <= k_pre; k++) {
142:     PetscCall(positiveVectorUpdate(rand, x, f));
143:     PetscCall(MatLMVMUpdate(A, x, f));
144:   }
145:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_copy));
146:   PetscCall(testMatEqual(rand, A, A_copy, &equal));
147:   PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatCopy() not the same after initial copy");

149:   PetscCall(VecDuplicate(x, &x_copy));
150:   PetscCall(VecCopy(x, x_copy));
151:   PetscCall(VecDuplicate(f, &f_copy));
152:   PetscCall(VecCopy(f, f_copy));

154:   PetscCall(testUnchangedBegin(rand, A_copy, &v1, &v2, &v3));
155:   PetscCall(positiveVectorUpdate(rand, x, f));
156:   PetscCall(MatLMVMUpdate(A, x, f));
157:   PetscCall(testUnchangedEnd(rand, A_copy, &v1, &v2, &v3, &equal));
158:   PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatLMVMUpdate() to original matrix affects copy");

160:   PetscCall(testUnchangedBegin(rand, A, &v1, &v2, &v3));
161:   PetscCall(positiveVectorUpdate(rand, x_copy, f_copy));
162:   PetscCall(MatLMVMUpdate(A_copy, x_copy, f_copy));
163:   PetscCall(testUnchangedEnd(rand, A, &v1, &v2, &v3, &equal));
164:   PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatLMVMUpdate() to copy matrix affects original");

166:   PetscCall(VecDestroy(&f_copy));
167:   PetscCall(VecDestroy(&x_copy));
168:   PetscCall(MatDestroy(&A_copy));
169:   PetscCall(MatDestroy(&A));
170:   PetscCall(VecDestroy(&f));
171:   PetscCall(VecDestroy(&x));
172:   PetscCall(VecDestroy(&u));
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: int main(int argc, char **argv)
177: {
178:   MPI_Comm    comm;
179:   PetscRandom rand;

181:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
182:   comm = PETSC_COMM_WORLD;
183:   PetscCall(PetscRandomCreate(comm, &rand));
184:   PetscCall(PetscRandomSetFromOptions(rand));
185:   PetscCall(KSPInitializePackage());
186:   PetscCall(testMatLMVMCopy(rand));
187:   PetscCall(PetscRandomDestroy(&rand));
188:   PetscCall(PetscFinalize());
189:   return 0;
190: }

192: /*TEST

194:   test:
195:     suffix: 0
196:     args: -mat_type {{lmvmbfgs lmvmdfp lmvmsr1 lmvmbroyden lmvmbadbroyden lmvmsymbroyden lmvmsymbadbroyden lmvmdiagbroyden lmvmdbfgs lmvmddfp lmvmdqn}}

198: TEST*/