aziis98
/
arnoldi-distribuito
mirror of https://github.com/aziis98/arnoldi-distribuito
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

feat: this maybe works

Browse Source
main
Antonio De Lucreziis 9 months ago
parent 506acf04e3
commit a1937f7edd
4 changed files with 116 additions and 60 deletions
Show Stats Download Patch File Download Diff File

3
.gitignore vendored
Unescape Escape View File

@ -4,3 +4,6 @@
# Binaries
build/
# Clang
.cache/

160
main.c
Unescape Escape View File

@ -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;

13
matrices/laplacian/julia-laplaciano.jl
Unescape Escape View File

@ -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))

BIN
matrices/laplacian/laplacian-discretization-3d.mat
View File

Binary file not shown.
Write Preview
Loading…
Cancel
Save