1 const char help[] = "Test that MatCopy() does not affect the copied LMVM matrix"; 2 3 #include <petsc.h> 4 5 static PetscErrorCode positiveVectorUpdate(PetscRandom rand, Vec x, Vec f) 6 { 7 Vec _x, _f; 8 PetscScalar dot; 9 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 } 22 23 static PetscErrorCode VecEqualToTolerance(Vec a, Vec b, NormType norm_type, PetscReal tol, PetscBool *flg) 24 { 25 Vec diff; 26 PetscReal diff_norm; 27 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 } 37 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; 42 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; 57 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; 64 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 } 74 75 static PetscErrorCode testUnchangedBegin(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z) 76 { 77 Vec _x, _y, _z; 78 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 } 90 91 static PetscErrorCode testUnchangedEnd(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z, PetscBool *unchanged) 92 { 93 Vec _x, _y, _z, _y2, _z2; 94 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 } 116 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; 126 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"); 148 149 PetscCall(VecDuplicate(x, &x_copy)); 150 PetscCall(VecCopy(x, x_copy)); 151 PetscCall(VecDuplicate(f, &f_copy)); 152 PetscCall(VecCopy(f, f_copy)); 153 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"); 159 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"); 165 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 } 175 176 int main(int argc, char **argv) 177 { 178 MPI_Comm comm; 179 PetscRandom rand; 180 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 } 191 192 /*TEST 193 194 test: 195 suffix: 0 196 args: -mat_type {{lmvmbfgs lmvmdfp lmvmsr1 lmvmbroyden lmvmbadbroyden lmvmsymbroyden lmvmsymbadbroyden lmvmdiagbroyden lmvmdbfgs lmvmddfp lmvmdqn}} 197 198 TEST*/ 199