From f0444e23992f96400a2a65b2131612d056f9afe7 Mon Sep 17 00:00:00 2001 From: Antonio De Lucreziis Date: Sat, 22 Mar 2025 03:50:59 +0100 Subject: [PATCH] fixed final hesseberg bug --- .gitignore | 4 +- arnoldi.c | 84 +++++++++++++++++++++++++++------ matrices/laplacian/eigs.jl | 20 ++++---- matrices/laplacian/laplacian.jl | 2 +- 4 files changed, 85 insertions(+), 25 deletions(-) diff --git a/.gitignore b/.gitignore index 2515d99..7d7c26d 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,6 @@ build*/ *.log *.out -*.lock? \ No newline at end of file +*.lock? + +*.mat \ No newline at end of file diff --git a/arnoldi.c b/arnoldi.c index a811d92..e6ad12b 100644 --- a/arnoldi.c +++ b/arnoldi.c @@ -1,3 +1,5 @@ +#include + #include #include @@ -17,6 +19,7 @@ #include #include #include +#include #include static char help[] = "Example PETSc program\n\n"; @@ -129,7 +132,20 @@ int main(int argc, char **argv) { VecSetSizes(b, PETSC_DECIDE, n); VecSetType(b, VECMPI); - VecSet(b, 1.0); + // VecSet(b, 1.0); + + // seed random number generator + srand((unsigned int)time(NULL)); + + // fill b with random values + for (PetscInt i = 0; i < n; i++) { + double val = (double)rand() / RAND_MAX; + VecSetValue(b, i, val, INSERT_VALUES); + } + + VecAssemblyBegin(b); + VecAssemblyEnd(b); + // VecSetValue(b, 0, 1.0, INSERT_VALUES); // MatView(A, PETSC_VIEWER_STDOUT_WORLD); @@ -156,14 +172,18 @@ int main(int argc, char **argv) { // MatSetType(H, MATMPIAIJ); - double *H = (double *)malloc((l + 1) * l * sizeof(double)); + double *h = (double *)malloc((l + 1) * l * sizeof(double)); + + for (PetscInt ii = 0; ii < (l + 1) * l; ii++) { + h[ii] = 0.0; + } // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Starting iteration\n")); PetscLogDouble arnoldi_start_time; PetscCall(PetscTime(&arnoldi_start_time)); { - PetscCall(ArnoldiIteration(A, b, l, n, Q, H)); + PetscCall(ArnoldiIteration(A, b, l, n, Q, h)); } PetscLogDouble arnoldi_end_time; PetscCall(PetscTime(&arnoldi_end_time)); @@ -181,6 +201,14 @@ int main(int argc, char **argv) { // PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 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("%.2f\t", h[i * (l + 1) + j]); + } + printf("\n"); + } double *wr = (double *)malloc(l * sizeof(double)); double *wi = (double *)malloc(l * sizeof(double)); @@ -192,8 +220,9 @@ int main(int argc, char **argv) { PetscCall(PetscTime(&lapack_start_time)); { // call LAPACK function "DHSEQR" to compute the eigenvalues of the Hessenberg matrix - LAPACKE_dhseqr(LAPACK_ROW_MAJOR, 'E', 'I', l, 1, l, H, l, wr, wi, z, l); + LAPACKE_dhseqr(LAPACK_ROW_MAJOR, 'E', 'I', l, 1, l, h, l + 1, wr, wi, z, l); } + PetscLogDouble lapack_end_time; PetscCall(PetscTime(&lapack_end_time)); PetscLogDouble lapack_time = lapack_end_time - lapack_start_time; @@ -203,7 +232,7 @@ int main(int argc, char **argv) { printf("H = \n"); for (int i = 0; i < l + 1; i++) { for (int j = 0; j < l; j++) { - printf("%.2f ", H[i * (l + 1) + j]); + printf("%.2f\t", h[i * (l + 1) + j]); } printf("\n"); } @@ -254,11 +283,29 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, do Vec q; - for (PetscInt i = 0; i < n + 1; i++) { - for (PetscInt j = 0; j < n; j++) { - h[i * (n + 1) + j] = 0.0; - } - } + if (rank == 0) + printf("n = %d, m = %d\n", n, m); + + /* + + eps = 1e-12 + h = np.zeros((n + 1, n)) + Q = np.zeros((A.shape[0], n + 1)) + # Normalize the input vector + Q[:, 0] = b / np.linalg.norm(b, 2) # Use it as the first Krylov vector + for k in range(1, n + 1): + v = np.dot(A, Q[:, k - 1]) # Generate a new candidate vector + for j in range(k): # Subtract the projections on previous vectors + h[j, k - 1] = np.dot(Q[:, j].conj(), v) + v = v - h[j, k - 1] * Q[:, j] + h[k, k - 1] = np.linalg.norm(v, 2) + if h[k, k - 1] > eps: # Add the produced vector to the list, unless + Q[:, k] = v / h[k, k - 1] + else: # If that happens, stop iterating. + return Q, h + + + */ PetscCall(VecDuplicate(b, &q)); PetscCall(VecCopy(b, q)); @@ -278,15 +325,20 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, do PetscCall(MatMult(A, Q[k - 1], v)); // Reorthogonalization using modified Gram-Schmidt - for (PetscInt j = 0; j < k - 1; j++) { // anche solo 3 + // fill all rows from 0 to k - 1 at column k - 1 + for (PetscInt j = 0; j < k; j++) { // anche solo 3 + PetscScalar h_ij; // h_(j, k - 1) = q_j . v PetscCall(VecDot(Q[j], v, &h_ij)); - + + if (rank == 0) + printf("h[%d, %d] = %f\n", j, k - 1, h_ij); + + h[j * (n + 1) + k - 1] = h_ij; + // v -= h_ij * q_j PetscCall(VecAXPY(v, -h_ij, Q[j])); - - h[j * (n + 1) + k - 1] = h_ij; } // Normalize: @@ -296,6 +348,10 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, do PetscScalar h_ij; PetscCall(VecNorm(v, NORM_2, &h_ij)); + + if (rank == 0) + printf("h[%d, %d] = %f\n", k, k - 1, h_ij); + h[k * (n + 1) + k - 1] = h_ij; // Check for convergence diff --git a/matrices/laplacian/eigs.jl b/matrices/laplacian/eigs.jl index 73cd409..4d83012 100644 --- a/matrices/laplacian/eigs.jl +++ b/matrices/laplacian/eigs.jl @@ -11,15 +11,17 @@ ex = fill(1, nx) ey = fill(1, ny) ez = fill(1, nz) -Dxx = spdiagm(-1 => ex, 0 => -2 * ex, +1 => ex) -Dyy = spdiagm(-1 => ey, 0 => -2 * ey, +1 => ey) -Dzz = spdiagm(-1 => ez, 0 => -2 * ez, +1 => ez) +Dxx = diagm(-1 => ex, 0 => -2 * ex, +1 => ex) +Dyy = diagm(-1 => ey, 0 => -2 * ey, +1 => ey) +Dzz = diagm(-1 => ez, 0 => -2 * ez, +1 => ez) -Ix = spdiagm(0 => [ex; 1]) -Iy = spdiagm(0 => [ey; 1]) -Iz = spdiagm(0 => [ez; 1]) +Ix = diagm(0 => [ex; 1]) +Iy = diagm(0 => [ey; 1]) +Iz = diagm(0 => [ez; 1]) -# L = kron(Dxx, Iy, Iz) + kron(Ix, Dyy, Iz) + kron(Ix, Iy, Dzz) +L = kron(Dxx, Iy, Iz) + kron(Ix, Dyy, Iz) + kron(Ix, Iy, Dzz) + +# display(eigvals(L)) # # 10 x 17 grid # nx = 11 - 1 @@ -40,7 +42,7 @@ println("Laplacian matrix L:") display(sparse(L)) -l = 100 +l = 400 # arnoldi iteration Q = Matrix{Float64}(undef, size(L, 1), l) @@ -66,6 +68,6 @@ end H = H[1:(l-1), 1:(l-1)] println("Hessenberg matrix H:") -display(H) +display(sparse(H)) display(eigvals(H)) diff --git a/matrices/laplacian/laplacian.jl b/matrices/laplacian/laplacian.jl index 3de2fd9..e7fe4f3 100644 --- a/matrices/laplacian/laplacian.jl +++ b/matrices/laplacian/laplacian.jl @@ -1,7 +1,7 @@ using SparseArrays using MAT -if len(ARGS) == 0 +if length(ARGS) == 0 println("Usage: julia laplacian.jl ") exit(1) end