Solve linear equations A * X = B
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=double), | intent(in), | dimension(:,:) | :: | a | ||
| real(kind=double), | intent(in), | dimension(:) | :: | b | ||
| real(kind=double), | intent(out), | dimension(:) | :: | x | 
  subroutine linsolve(a,b,x)
   !+ Solve linear equations A * X = B
   real(double), dimension(:,:), intent(IN)  :: a ! assume a square
   real(double), dimension(:), intent(IN)    :: b ! assume size(b) == size(a,1)
   real(double), dimension(:), intent(OUT)   :: x ! assume size(x) == size(b)
#ifdef _USE_BLAS
   real(double), dimension(size(b),size(b)) :: lu
   integer, dimension(size(b)) :: ipiv
   integer :: n,info
   n = size(b)
   ! first factorize a
   lu(:,:) = a(:,:)
   call DGETRF(n,n,lu,n,ipiv,info)
   if (info /= 0) stop 'sub linsolve: DGETRF failed!'
   x(:) = b(:)
   call DGETRS('N',n,1,lu,n,ipiv,x,n,info)
   if (info /= 0) stop 'sub linsolve: DGETRS failed!'
#else
   real(double), dimension(size(b),size(b)) :: a_inv
   call matinv(a, a_inv)
   x = matmul(a_inv, b)!coeffs determined from states at t=0
#endif
  end subroutine linsolve