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

finito?

Browse Source
next
Antonio De Lucreziis 1 year ago
parent 4f1adee70d
commit ee8054b7f8
6 changed files with 244 additions and 97 deletions
Show Stats Download Patch File Download Diff File

1
.clang-format
Unescape Escape View File

@ -1,2 +1,3 @@
# .clang-format
IndentWidth: 4
ColumnLimit: 120

4
CMakeLists.txt
Unescape Escape View File

@ -24,3 +24,7 @@ target_link_libraries(main PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConf
add_executable(arnoldi arnoldi.c)
target_include_directories(arnoldi PUBLIC ${PETSc_INCLUDE_DIRS} ${OpenBLAS_INCLUDE_DIRS} ${OpenMPI_INCLUDE_DIRS})
target_link_libraries(arnoldi PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConfig::OpenMPI)
add_executable(ex10 ex10.c)
target_include_directories(ex10 PUBLIC ${PETSc_INCLUDE_DIRS} ${OpenBLAS_INCLUDE_DIRS} ${OpenMPI_INCLUDE_DIRS})
target_link_libraries(ex10 PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConfig::OpenMPI)

236
arnoldi.c
Unescape Escape View File

@ -1,20 +1,41 @@
#include <lapack.h>
#include <lapacke.h>
#include <mpi.h>
#include <petscerror.h>
#include <petscmat.h>
#include <petscsys.h>
#include <petscsystypes.h>
#include <petscvec.h>
#include <petscviewer.h>
#include <petscviewerhdf5.h>
#include <petscdm.h>
#include <petscdmda.h>
#include <petscksp.h>
#include <petscviewertypes.h>
#include <stdlib.h>
static char help[] = "Example PETSc program\n\n";
// extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
// extern PetscErrorCode ComputeRHS(KSP, Vec, void *);
// extern PetscErrorCode ComputeInitialSolution(DM, Vec);
void swap(double *a, double *b) {
double t = *a;
*a = *b;
*b = t;
}
// int compare_double(const void *a, const void *b) {
// int aa = *(int *)a;
// int bb = *(int *)b;
PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H);
// if (aa < bb)
// return -1;
// if (aa > bb)
// return 1;
// return 0;
// }
PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, double *h);
int main(int argc, char **argv) {
Vec b;
@ -26,78 +47,47 @@ int main(int argc, char **argv) {
PetscBool flg;
PetscOptionsGetInt(NULL, NULL, "-n", &n, &flg);
if (!flg)
n = 176;
n = -1;
PetscOptionsGetInt(NULL, NULL, "-l", &l, &flg);
if (!flg)
l = 4;
VecCreate(PETSC_COMM_WORLD, &b);
VecSetSizes(b, PETSC_DECIDE, n);
VecSetType(b, VECMPI);
VecSet(b, 1.0);
// VecSetValue(b, 0, 1.0, INSERT_VALUES);
l = 10;
Mat A;
MatCreate(PETSC_COMM_WORLD, &A);
// 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(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);
VecView(b, PETSC_VIEWER_DRAW_WORLD);
VecView(b, PETSC_VIEWER_STDOUT_WORLD);
// TIME: Assembly construction
printf("[Arnoldi] Allocating memory for Krylov subspace basis\n");
if (n == -1) {
MatGetSize(A, &n, NULL);
}
VecCreate(PETSC_COMM_WORLD, &b);
VecSetSizes(b, PETSC_DECIDE, n);
VecSetType(b, VECMPI);
VecSet(b, 1.0);
// VecSetValue(b, 0, 1.0, INSERT_VALUES);
// MatView(A, PETSC_VIEWER_STDOUT_WORLD);
// VecView(b, PETSC_VIEWER_STDOUT_WORLD);
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Allocating memory for Krylov subspace basis\n"));
Vec *Q;
PetscMalloc1(l, &Q);
@ -106,65 +96,110 @@ int main(int argc, char **argv) {
VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]);
}
printf("[Arnoldi] Constructing Hessenberg matrix\n");
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Constructing Hessenberg matrix\n"));
// Mat H;
// PetscCall(MatCreate(PETSC_COMM_SELF, &H));
// PetscCall(PetscObjectSetName((PetscObject)H, "hessenberg"));
// PetscCall(MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l));
// PetscCall(MatSetType(H, MATDENSE));
// PetscCall(MatCreate(PETSC_COMM_WORLD, &H));
Mat H;
PetscCall(MatCreate(PETSC_COMM_SELF, &H));
PetscCall(MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l));
PetscCall(MatSetType(H, MATDENSE));
// MatSetType(H, MATMPIAIJ);
printf("[Arnoldi] Starting iteration\n");
double *H = (double *)malloc((l + 1) * l * sizeof(double));
PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Starting iteration\n"));
PetscCall(ArnoldiIteration(A, b, l, Q, H));
// TIME:START
PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
PetscCall(ArnoldiIteration(A, b, l, n, Q, H));
printf("[Arnoldi] Done\n");
// PetscCall(MatSetValue(H, 2, 3, -1, INSERT_VALUES));
// PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
// PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
MatView(H, PETSC_VIEWER_DRAW_WORLD);
MatView(H, PETSC_VIEWER_STDOUT_WORLD);
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Done\n"));
// print Hessenberg matrix to file
int rank;
PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
PetscViewer v2;
PetscCall(PetscViewerCreate(PETSC_COMM_SELF, &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));
if (rank == 0) {
// print Hessenberg matrix
printf("H = \n");
for (int i = 0; i < l + 1; i++) {
for (int j = 0; j < l; j++) {
printf("%.3f ", H[i * l + j]);
}
printf("\n");
}
// for (PetscInt i = 0; i < l + 1; i++) {
// VecView(Q[i], PETSC_VIEWER_STDOUT_WORLD);
// }
// call LAPACK function "DHSEQR" to compute the eigenvalues of the Hessenberg matrix
double *wr = (double *)malloc(l * sizeof(double));
double *wi = (double *)malloc(l * sizeof(double));
double *z = (double *)malloc(l * l * sizeof(double));
double *work = (double *)malloc(3 * l * sizeof(double));
int info;
LAPACKE_dhseqr(LAPACK_ROW_MAJOR, 'E', 'I', l, 1, l, H, l, wr, wi, z, l);
// TIME:END
// sort eigenvalues
for (int i = 0; i < l; i++) {
for (int j = i + 1; j < l; j++) {
if (wr[i] > wr[j]) {
swap(&wr[i], &wr[j]);
swap(&wi[i], &wi[j]);
for (int k = 0; k < l; k++) {
swap(&z[i * l + k], &z[j * l + k]);
}
}
}
}
// MatView(H, PETSC_VIEWER_STDOUT_WORLD);
// print eigenvalues
printf("Eigenvalues = \n");
for (int i = 0; i < l; i++) {
printf("%.3f + %.3f i\n", wr[i], wi[i]);
}
}
// for (PetscInt i = 0; i < l + 1; i++) {
// PetscCall(VecDestroy(&Q[i]));
// print eigenvectors
// printf("Eigenvectors = \n");
// for (int i = 0; i < l; i++) {
// for (int j = 0; j < l; j++) {
// printf("%f ", z[i * l + j]);
// }
// printf("\n");
// }
// PetscCall(PetscFree(Q));
PetscCall(MatDestroy(&A));
PetscCall(VecDestroy(&b));
// PetscCall(MatDestroy(&H));
PetscFinalize();
return 0;
}
PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) {
PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, double *h) {
PetscFunctionBeginUser;
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
printf("Rank %d\n", rank);
PetscScalar eps = 1e-12;
PetscInt m;
PetscCall(VecGetSize(b, &m));
Vec q;
PetscCall(MatZeroEntries(H));
for (PetscInt i = 0; i < n + 1; i++) {
for (PetscInt j = 0; j < n; j++) {
h[i * n + j] = 0.0;
}
}
PetscCall(VecDuplicate(b, &q));
PetscCall(VecCopy(b, q));
@ -173,32 +208,39 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) {
Q[0] = q;
for (PetscInt k = 1; k < n + 1; k++) {
// printf("[Arnoldi] Iteration %d\n", k);
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Iteration %d\n", k));
Vec 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);
for (PetscInt j = 0; j < k; j++) { // anche solo 3
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Reorthogonalization %d %d\n", k, j));
PetscScalar h_ij;
PetscCall(VecDot(Q[j], v, &h_ij));
PetscScalar h;
PetscCall(VecDot(Q[j], v, &h));
PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES));
PetscCall(VecAXPY(v, -h, Q[j]));
// PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES));
h[j * n + k - 1] = h_ij;
PetscCall(VecAXPY(v, -h_ij, Q[j]));
}
// Normalize
PetscScalar h;
PetscCall(VecNorm(v, NORM_2, &h));
PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES));
PetscScalar h_ij;
PetscCall(VecNorm(v, NORM_2, &h_ij));
// PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES));
h[k * n + k - 1] = h_ij;
// Check for convergence
if (h > eps) {
if (h_ij > eps) {
PetscCall(VecNormalize(v, NULL));
Q[k] = v;
} else {
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Early breakdown"));
break;
}
}

82
ex10.c
Unescape Escape View File

@ -0,0 +1,82 @@
#include <petscdm.h>
#include <petscdmda.h>
#include <petscsys.h>
#include <petscviewerhdf5.h>
static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n";
int main(int argc, char **argv) {
DM da2D;
PetscInt i, j, ixs, ixm, iys, iym;
PetscViewer H5viewer;
PetscScalar xm = -1.0, xp = 1.0;
PetscScalar ym = -1.0, yp = 1.0;
PetscScalar value = 1.0, dx, dy;
PetscInt Nx = 40, Ny = 40;
Vec gauss, input;
PetscScalar **gauss_ptr;
PetscReal norm;
const char *vecname;
dx = (xp - xm) / (Nx - 1);
dy = (yp - ym) / (Ny - 1);
/* Initialize the Petsc context */
PetscFunctionBeginUser;
PetscCall(PetscInitialize(&argc, &argv, NULL, help));
PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, Nx, Ny,
PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da2D));
PetscCall(DMSetFromOptions(da2D));
PetscCall(DMSetUp(da2D));
/* Set the coordinates */
PetscCall(DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0));
/* Declare gauss as a DMDA component */
PetscCall(DMCreateGlobalVector(da2D, &gauss));
PetscCall(PetscObjectSetName((PetscObject)gauss, "pressure"));
/* Initialize vector gauss with a constant value (=1) */
PetscCall(VecSet(gauss, value));
/* Get the coordinates of the corners for each process */
PetscCall(DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0));
/* Build the gaussian profile (exp(-x^2-y^2)) */
PetscCall(DMDAVecGetArray(da2D, gauss, &gauss_ptr));
for (j = iys; j < iys + iym; j++) {
for (i = ixs; i < ixs + ixm; i++)
gauss_ptr[j][i] = PetscExpScalar(-(xm + i * dx) * (xm + i * dx) - (ym + j * dy) * (ym + j * dy));
}
PetscCall(DMDAVecRestoreArray(da2D, gauss, &gauss_ptr));
/* Create the HDF5 viewer */
PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer));
PetscCall(PetscViewerSetFromOptions(H5viewer));
/* Write the H5 file */
PetscCall(VecView(gauss, H5viewer));
/* Close the viewer */
PetscCall(PetscViewerDestroy(&H5viewer));
PetscCall(VecDuplicate(gauss, &input));
PetscCall(PetscObjectGetName((PetscObject)gauss, &vecname));
PetscCall(PetscObjectSetName((PetscObject)input, vecname));
/* Create the HDF5 viewer for reading */
PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer));
PetscCall(PetscViewerSetFromOptions(H5viewer));
PetscCall(VecLoad(input, H5viewer));
PetscCall(PetscViewerDestroy(&H5viewer));
PetscCall(VecAXPY(input, -1.0, gauss));
PetscCall(VecNorm(input, NORM_2, &norm));
PetscCheck(norm <= 1.e-6, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Vec read in does not match vector written out");
PetscCall(VecDestroy(&input));
PetscCall(VecDestroy(&gauss));
PetscCall(DMDestroy(&da2D));
PetscCall(PetscFinalize());
return 0;
}

1
main.c
Unescape Escape View File

@ -107,6 +107,7 @@ int main(int argc, char **argv) {
PetscCall(
PetscPrintf(PETSC_COMM_WORLD, "Vector length %" PetscInt_FMT "\n", n));
PetscCall(VecMax(x, &maxind, &maxval));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,
"VecMax %g, VecInd %" PetscInt_FMT "\n",

17
matrices/laplacian/eigs.jl
Unescape Escape View File

@ -0,0 +1,17 @@
using SparseArrays
using LinearAlgebra
using MAT
# 11 x 16
nx = 10
ny = 15
ex = fill(1, nx)
ey = fill(1, ny)
Dxx = diagm(-1 => ex, 0 => -2 * ex, +1 => ex)
Dyy = diagm(-1 => ey, 0 => -2 * ey, +1 => ey)
L = kron(Dyy, diagm(0 => [ex; 1])) + kron(diagm(0 => [ey; 1]), Dxx);
println("Laplacian matrix L:")
display(sparse(L))
display(eigvals(L))
Write Preview
Loading…
Cancel
Save