diff --git a/.gitignore b/.gitignore index 57fd4da..08d45fa 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,6 @@ # Binaries build/ + +# Clang +.cache/ \ No newline at end of file diff --git a/main.c b/main.c index 9fca88c..a17041e 100644 --- a/main.c +++ b/main.c @@ -17,7 +17,6 @@ static char help[] = "Example PETSc program\n\n"; PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H); int main(int argc, char **argv) { - Mat A; Vec b; PetscInt n, l; @@ -27,7 +26,7 @@ int main(int argc, char **argv) { PetscBool flg; PetscOptionsGetInt(NULL, NULL, "-n", &n, &flg); if (!flg) - n = 10; + n = 176; PetscOptionsGetInt(NULL, NULL, "-l", &l, &flg); if (!flg) @@ -40,42 +39,65 @@ int main(int argc, char **argv) { VecSet(b, 1.0); // VecSetValue(b, 0, 1.0, INSERT_VALUES); + Mat A; MatCreate(PETSC_COMM_WORLD, &A); - MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n); - - MatSetType(A, MATMPIAIJ); - - // A := diag(-1, 2, -1) - for (PetscInt i = 0; i < n; i++) { - PetscScalar v[3] = {-1.0, 2.0, -1.0}; - PetscInt col[3] = {i - 1, i, i + 1}; - PetscInt ncol = 0; - if (i > 0) { - col[ncol] = i - 1; - v[ncol] = -1.0; - ncol++; - } - col[ncol] = i; - v[ncol] = 2.0; - ncol++; - if (i < n - 1) { - col[ncol] = i + 1; - v[ncol] = -1.0; - ncol++; - } - MatSetValues(A, 1, &i, ncol, col, v, INSERT_VALUES); - // MatSetValue(A, i, i, (PetscScalar)(i + 1), INSERT_VALUES); - } + // MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n); + // MatSetType(A, MATMPIAIJ); + PetscViewer v; + PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &v)); + PetscCall(PetscViewerSetType(v, PETSCVIEWERHDF5)); + PetscCall(PetscViewerPushFormat(v, PETSC_VIEWER_HDF5_MAT)); + PetscCall(PetscViewerSetFromOptions(v)); + PetscCall(PetscViewerFileSetMode(v, FILE_MODE_READ)); + PetscCall(PetscViewerFileSetName( + v, "../matrices/laplacian/laplacian-discretization-3d.mat")); + + PetscCall(MatSetOptionsPrefix(A, "a_")); + PetscCall(PetscObjectSetName((PetscObject)A, "A")); + + // PetscCall( + // PetscOptionsGetString(NULL, NULL, "-f", name, sizeof(name), &flg)); + // PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, + // "Must provide a binary file for the matrix"); + + // // PetscCall(MatLoad(A, v)); + // PetscCall(PetscViewerBinaryOpen( + // PETSC_COMM_WORLD, + // "../matrices/laplacian/laplacian-discretization-3d.mat", + // FILE_MODE_READ, &v)); + PetscCall(MatLoad(A, v)); + + // // A := diag(-1, 2, -1) + // for (PetscInt i = 0; i < n; i++) { + // PetscScalar v[3] = {-1.0, 2.0, -1.0}; + // PetscInt col[3] = {i - 1, i, i + 1}; + // PetscInt ncol = 0; + // if (i > 0) { + // col[ncol] = i - 1; + // v[ncol] = -1.0; + // ncol++; + // } + // col[ncol] = i; + // v[ncol] = 2.0; + // ncol++; + // if (i < n - 1) { + // col[ncol] = i + 1; + // v[ncol] = -1.0; + // ncol++; + // } + // MatSetValues(A, 1, &i, ncol, col, v, INSERT_VALUES); + // // MatSetValue(A, i, i, (PetscScalar)(i + 1), INSERT_VALUES); + // } MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); // MatView(A, PETSC_VIEWER_DRAW_WORLD); - MatView(A, PETSC_VIEWER_STDOUT_WORLD); + // MatView(A, PETSC_VIEWER_STDOUT_WORLD); // VecView(b, PETSC_VIEWER_DRAW_WORLD); - VecView(b, PETSC_VIEWER_STDOUT_WORLD); + // VecView(b, PETSC_VIEWER_STDOUT_WORLD); - printf("Allocating memory for Krylov subspace basis\n"); + printf("[Arnoldi] Allocating memory for Krylov subspace basis\n"); Vec *Q; PetscMalloc1(l, &Q); @@ -84,34 +106,48 @@ int main(int argc, char **argv) { VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]); } - printf("Constructing Hessenberg matrix\n"); + printf("[Arnoldi] Constructing Hessenberg matrix\n"); Mat H; - MatCreate(PETSC_COMM_WORLD, &H); - MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l); + PetscCall(MatCreate(PETSC_COMM_WORLD, &H)); + PetscCall(MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l)); + PetscCall(MatSetType(H, MATDENSE)); // MatSetType(H, MATMPIAIJ); - MatSetType(H, MATDENSE); - printf("Starting Arnoldi iteration\n"); + printf("[Arnoldi] Starting iteration\n"); + + PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); - MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY); PetscCall(ArnoldiIteration(A, b, l, Q, H)); - MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY); - for (PetscInt i = 0; i < l + 1; i++) { - VecView(Q[i], PETSC_VIEWER_STDOUT_WORLD); - } + PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); - MatView(H, PETSC_VIEWER_STDOUT_WORLD); + printf("[Arnoldi] Done\n"); - for (PetscInt i = 0; i < l + 1; i++) { - VecDestroy(&Q[i]); - } + // print Hessenberg matrix to file + + PetscViewer v2; + PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &v2)); + PetscCall(PetscViewerSetType(v2, PETSCVIEWERHDF5)); + PetscCall(PetscViewerPushFormat(v2, PETSC_VIEWER_HDF5_MAT)); + PetscCall(PetscViewerFileSetMode(v2, FILE_MODE_WRITE)); + PetscCall(PetscViewerFileSetName(v2, "hessenberg.mat")); + PetscCall(MatView(H, v2)); + + // for (PetscInt i = 0; i < l + 1; i++) { + // VecView(Q[i], PETSC_VIEWER_STDOUT_WORLD); + // } + + // MatView(H, PETSC_VIEWER_STDOUT_WORLD); - // PetscFree(Q); + // for (PetscInt i = 0; i < l + 1; i++) { + // PetscCall(VecDestroy(&Q[i])); + // } - MatDestroy(&A); - VecDestroy(&b); + // PetscCall(PetscFree(Q)); + + PetscCall(MatDestroy(&A)); + PetscCall(VecDestroy(&b)); PetscFinalize(); return 0; } @@ -121,39 +157,43 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) { PetscScalar eps = 1e-12; PetscInt m; - VecGetSize(b, &m); + PetscCall(VecGetSize(b, &m)); Vec q; - MatZeroEntries(H); + PetscCall(MatZeroEntries(H)); - VecDuplicate(b, &q); - VecCopy(b, q); - VecNormalize(q, NULL); + PetscCall(VecDuplicate(b, &q)); + PetscCall(VecCopy(b, q)); + PetscCall(VecNormalize(q, NULL)); Q[0] = q; for (PetscInt k = 1; k < n + 1; k++) { + // printf("[Arnoldi] Iteration %d\n", k); + Vec v; - VecDuplicate(b, &v); - MatMult(A, Q[k - 1], v); + PetscCall(VecDuplicate(b, &v)); + PetscCall(MatMult(A, Q[k - 1], v)); // Reorthogonalization using modified Gram-Schmidt for (PetscInt j = 0; j < k; j++) { + // printf("[Arnoldi] Reorthogonalization %d\n", j); + PetscScalar h; - VecDot(Q[j], v, &h); - MatSetValue(H, j, k - 1, h, INSERT_VALUES); - VecAXPY(v, -h, Q[j]); + PetscCall(VecDot(Q[j], v, &h)); + PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES)); + PetscCall(VecAXPY(v, -h, Q[j])); } // Normalize PetscScalar h; - VecNorm(v, NORM_2, &h); - MatSetValue(H, k, k - 1, h, INSERT_VALUES); + PetscCall(VecNorm(v, NORM_2, &h)); + PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES)); // Check for convergence if (h > eps) { - VecNormalize(v, NULL); + PetscCall(VecNormalize(v, NULL)); Q[k] = v; } else { break; diff --git a/matrices/laplacian/julia-laplaciano.jl b/matrices/laplacian/julia-laplaciano.jl new file mode 100644 index 0000000..dcc29bd --- /dev/null +++ b/matrices/laplacian/julia-laplaciano.jl @@ -0,0 +1,13 @@ +using SparseArrays +using MAT + +# 11 x 16 +nx = 10 +ny = 15 +ex = fill(1, nx) +ey = fill(1, ny) +Dxx = spdiagm(-1 => ex, 0 => -2 * ex, +1 => ex) +Dyy = spdiagm(-1 => ey, 0 => -2 * ey, +1 => ey) +L = kron(Dyy, spdiagm(0 => [ex; 1])) + kron(spdiagm(0 => [ey; 1]), Dxx); + +matwrite("laplacian-discretization-3d.mat", Dict("A" => L)) diff --git a/matrices/laplacian/laplacian-discretization-3d.mat b/matrices/laplacian/laplacian-discretization-3d.mat new file mode 100644 index 0000000..bad09f4 Binary files /dev/null and b/matrices/laplacian/laplacian-discretization-3d.mat differ