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 // unlike MatTestEqual(), this test tests MatMult() and MatSolve() 24 static PetscErrorCode testMatEqual(PetscRandom rand, Mat A, Mat B, PetscBool *flg) 25 { 26 Vec x, y_A, y_B; 27 28 PetscFunctionBegin; 29 *flg = PETSC_TRUE; 30 PetscCall(MatCreateVecs(A, &x, &y_A)); 31 PetscCall(MatCreateVecs(B, NULL, &y_B)); 32 PetscCall(VecSetRandom(x, rand)); 33 PetscCall(MatMult(A, x, y_A)); 34 PetscCall(MatMult(B, x, y_B)); 35 PetscCall(VecEqual(y_A, y_B, flg)); 36 if (*flg == PETSC_TRUE) { 37 PetscCall(MatSolve(A, x, y_A)); 38 PetscCall(MatSolve(B, x, y_B)); 39 PetscCall(VecEqual(y_A, y_B, flg)); 40 } 41 PetscCall(VecDestroy(&y_B)); 42 PetscCall(VecDestroy(&y_A)); 43 PetscCall(VecDestroy(&x)); 44 PetscFunctionReturn(PETSC_SUCCESS); 45 } 46 47 static PetscErrorCode testUnchangedBegin(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z) 48 { 49 Vec _x, _y, _z; 50 51 PetscFunctionBegin; 52 PetscCall(MatCreateVecs(A, &_x, &_y)); 53 PetscCall(MatCreateVecs(A, NULL, &_z)); 54 PetscCall(VecSetRandom(_x, rand)); 55 PetscCall(MatMult(A, _x, _y)); 56 PetscCall(MatSolve(A, _x, _z)); 57 *x = _x; 58 *y = _y; 59 *z = _z; 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 static PetscErrorCode testUnchangedEnd(PetscRandom rand, Mat A, Vec *x, Vec *y, Vec *z, PetscBool *unchanged) 64 { 65 Vec _x, _y, _z, _y2, _z2; 66 67 PetscFunctionBegin; 68 *unchanged = PETSC_TRUE; 69 _x = *x; 70 _y = *y; 71 _z = *z; 72 *x = NULL; 73 *y = NULL; 74 *z = NULL; 75 PetscCall(MatCreateVecs(A, NULL, &_y2)); 76 PetscCall(MatCreateVecs(A, NULL, &_z2)); 77 PetscCall(MatMult(A, _x, _y2)); 78 PetscCall(MatSolve(A, _x, _z2)); 79 PetscCall(VecEqual(_y, _y2, unchanged)); 80 if (*unchanged == PETSC_TRUE) PetscCall(VecEqual(_z, _z2, unchanged)); 81 PetscCall(VecDestroy(&_z2)); 82 PetscCall(VecDestroy(&_y2)); 83 PetscCall(VecDestroy(&_z)); 84 PetscCall(VecDestroy(&_y)); 85 PetscCall(VecDestroy(&_x)); 86 PetscFunctionReturn(PETSC_SUCCESS); 87 } 88 89 static PetscErrorCode testMatLMVMCopy(PetscRandom rand) 90 { 91 PetscInt N = 10; 92 MPI_Comm comm = PetscObjectComm((PetscObject)rand); 93 PetscInt k_pre = 2; // number of updates before copy 94 Mat A, A_copy; 95 Vec u, x, f, x_copy, f_copy, v1, v2, v3; 96 PetscBool equal; 97 PetscLayout layout; 98 99 PetscFunctionBegin; 100 PetscCall(VecCreateMPI(comm, PETSC_DECIDE, N, &u)); 101 PetscCall(VecSetFromOptions(u)); 102 PetscCall(VecSetUp(u)); 103 PetscCall(VecDuplicate(u, &x)); 104 PetscCall(VecDuplicate(u, &f)); 105 PetscCall(VecGetLayout(u, &layout)); 106 PetscCall(MatCreate(comm, &A)); 107 PetscCall(MatSetLayouts(A, layout, layout)); 108 PetscCall(MatSetType(A, MATLMVMBFGS)); 109 PetscCall(MatSetFromOptions(A)); 110 PetscCall(positiveVectorUpdate(rand, x, f)); 111 PetscCall(MatLMVMAllocate(A, x, f)); 112 PetscCall(MatSetUp(A)); 113 for (PetscInt k = 0; k <= k_pre; k++) { 114 PetscCall(positiveVectorUpdate(rand, x, f)); 115 PetscCall(MatLMVMUpdate(A, x, f)); 116 } 117 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_copy)); 118 PetscCall(testMatEqual(rand, A, A_copy, &equal)); 119 PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatCopy() not the same after initial copy"); 120 121 PetscCall(VecDuplicate(x, &x_copy)); 122 PetscCall(VecCopy(x, x_copy)); 123 PetscCall(VecDuplicate(f, &f_copy)); 124 PetscCall(VecCopy(f, f_copy)); 125 126 PetscCall(testUnchangedBegin(rand, A_copy, &v1, &v2, &v3)); 127 PetscCall(positiveVectorUpdate(rand, x, f)); 128 PetscCall(MatLMVMUpdate(A, x, f)); 129 PetscCall(testUnchangedEnd(rand, A_copy, &v1, &v2, &v3, &equal)); 130 PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatLMVMUpdate() to original matrix affects copy"); 131 132 PetscCall(testUnchangedBegin(rand, A, &v1, &v2, &v3)); 133 PetscCall(positiveVectorUpdate(rand, x_copy, f_copy)); 134 PetscCall(MatLMVMUpdate(A_copy, x_copy, f_copy)); 135 PetscCall(testUnchangedEnd(rand, A, &v1, &v2, &v3, &equal)); 136 PetscCheck(equal, comm, PETSC_ERR_PLIB, "MatLMVMUpdate() to copy matrix affects original"); 137 138 PetscCall(VecDestroy(&f_copy)); 139 PetscCall(VecDestroy(&x_copy)); 140 PetscCall(MatDestroy(&A_copy)); 141 PetscCall(MatDestroy(&A)); 142 PetscCall(VecDestroy(&f)); 143 PetscCall(VecDestroy(&x)); 144 PetscCall(VecDestroy(&u)); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 int main(int argc, char **argv) 149 { 150 MPI_Comm comm; 151 PetscRandom rand; 152 153 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 154 comm = PETSC_COMM_WORLD; 155 PetscCall(PetscRandomCreate(comm, &rand)); 156 PetscCall(PetscRandomSetFromOptions(rand)); 157 PetscCall(KSPInitializePackage()); 158 PetscCall(testMatLMVMCopy(rand)); 159 PetscCall(PetscRandomDestroy(&rand)); 160 PetscCall(PetscFinalize()); 161 return 0; 162 } 163 164 /*TEST 165 166 test: 167 suffix: 0 168 args: -mat_type {{lmvmbfgs lmvmdfp lmvmsr1 lmvmbroyden lmvmbadbroyden lmvmsymbroyden lmvmsymbadbroyden lmvmdiagbroyden lmvmdbfgs lmvmddfp lmvmdqn}} 169 170 TEST*/ 171