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 # dense != compact_dense 195 test: 196 suffix: 0 197 output_file: output/empty.out 198 args: -mat_lmvm_mult_algorithm {{recursive dense compact_dense}} -mat_type {{lmvmbfgs lmvmdfp lmvmbroyden lmvmbadbroyden}} 199 200 # dense == compact_dense 201 test: 202 suffix: 1 203 output_file: output/empty.out 204 args: -mat_lmvm_mult_algorithm {{recursive dense}} -mat_type {{lmvmsr1 lmvmsymbroyden lmvmsymbadbroyden}} 205 206 test: 207 suffix: 2 208 output_file: output/empty.out 209 args: -mat_type {{lmvmdiagbroyden lmvmdqn}} 210 211 test: 212 suffix: 3 213 output_file: output/empty.out 214 args: -mat_type lmvmdbfgs -mat_lbfgs_recursive {{0 1}} 215 216 test: 217 suffix: 4 218 output_file: output/empty.out 219 args: -mat_type lmvmddfp -mat_ldfp_recursive {{0 1}} 220 221 TEST*/ 222