disabling stdlib_linalg for now, since it takes too much to compiler.

app/critical.f90

@@ -2,7 +2,6 @@ program main
    !! Test the calculation of critical lines
    use yaeos
    use yaeos__math, only: solve_system
-   use stdlib_linalg, only: eigh
    use yaeos__equilibria_critical, only: &
       lambda1, F_critical, df_critical, CriticalLine, critical_line, critical_point
    implicit none

src/equilibria/critical.f90

@@ -197,7 +197,7 @@ real(pr) function lambda1(model, X, s, z0, zi, u, u_new, P)
       !! \[
       !!  \lambda_1(s) = \frac{d^2tpd}{ds^2} = 0
       !! \]
-      use stdlib_linalg, only: eigh, linalg_state_type
+      use yaeos__math_linalg, only: eigen
       class(ArModel), intent(in) :: model
       real(pr), intent(in) :: z0(:) !! Molar fractions of the first fluid
       real(pr), intent(in) :: zi(:) !! Molar fractions of the second fluid
@@ -216,7 +216,7 @@ real(pr) function lambda1(model, X, s, z0, zi, u, u_new, P)
       real(pr) :: dlnf_dn(size(z0), size(z0))
       real(pr) :: lambda(size(z0)), vectors(size(z0), size(z0))
 
-      type(linalg_state_type) :: stat
+      ! type(linalg_state_type) :: stat
 
       integer :: i, j, nc
       real(pr) :: M(size(z0), size(z0)), z(size(z0)), Pin
@@ -236,10 +236,7 @@ real(pr) function lambda1(model, X, s, z0, zi, u, u_new, P)
          end do
       end do
 
-      call eigh(A=M, lambda=lambda, vectors=vectors, err=stat)
-      if (.not. stat%ok()) then
-         write(*, *) stat%print_msg()
-      end if
+      call eigen(A=M, eigenvalues=lambda, eigenvectors=vectors)
 
       lambda1 = lambda(minloc(abs(lambda), dim=1))
       if (present(u_new)) u_new = vectors(:, minloc(abs(lambda), dim=1))
@@ -393,11 +390,11 @@ type(EquilibriumState) function critical_point(&
          df = df_critical(model, X, ns, S, z0, zi, u)
          dX = solve_system(A=df, b=-F)
 
-         do while(maxval(abs(dX/X)) > 1e-1)
+         do while(maxval(abs(dX)) > 1e-1)
             dX = dX/10
          end do
 
-         if (maxval(abs(X)) < 1e-5) exit
+         if (maxval(abs(F)) < 1e-5) exit
 
          X = X + dX
          l = lambda1(model, X, 0.0_pr, z0, zi, u, u_new)

src/math/linalg.f90

@@ -149,4 +149,53 @@ subroutine cubic_roots_rosendo(parameters, real_roots, complex_roots, flag)
       end if
    end subroutine
 
+   subroutine eigen(A, eigenvalues, eigenvectors)
+      !! # eigen
+      !!
+      !! # Description
+      !! Calculate the eigenvalues and eigenvectors of a real symmetric matrix
+      !! `A` using LAPACK's `dsyev`. The eigenvectors are stored in the columns
+      !! of `eigenvectors`. The eigenvalues are stored in `eigenvalues`.
+
+      ! use stdlib_linalg, only: eigh, linalg_state_type
+      ! type(linalg_state_type) :: stat
+      real(pr), intent(in out) :: A(:, :)
+      real(pr), intent(out) :: eigenvalues(:)
+      real(pr), optional, intent(out) :: eigenvectors(:, :)
+
+      interface
+         subroutine dsyev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO) 	
+            character(1) :: JOBZ
+            character(1) :: UPLO
+            integer :: N
+            double precision :: A(lda, *)
+            integer :: lda
+            double precision :: W(*)
+            double precision :: WORK(*)
+            integer :: LWORK
+            integer :: INFO
+         end subroutine
+      end interface
+      integer :: istat
+      real(pr) :: work(size(A, 1)*10)
+      real(pr) :: Ain(size(A, 1), size(A, 1))
+      integer :: n
+
+      n = size(A, 1)
+      Ain = A
+
+      call dsyev(jobz='V', uplo='U', n=n, A=Ain, lda=n, w=eigenvalues, &
+         work=work, lwork=size(work, 1), info=istat &
+      )
+      eigenvectors = Ain
+
+      ! call eigh(&
+      !    A=A, lambda=eigenvalues, vectors=eigenvectors, err=stat  &
+      ! )
+
+      ! if (.not. stat%ok()) then
+      !    write(*, *) stat%print_msg()
+      ! end if
+   end subroutine
+
 end module yaeos__math_linalg

test/test_critical.f90

@@ -83,6 +83,7 @@ program main
    end do
 
 contains
+
    type(CubicEoS) function get_model()
       real(pr) :: tc(nc), pc(nc), w(nc), kij(nc,nc)
       z0=[0.0656,0.3711,0.0538,0.0373,0.0261,0.0187,&