適当にソースコードを拾ってくる．

wget http://www.netlib.org/templates/double/BiCGSTAB.f
wget http://www.netlib.org/templates/double/BiCGSTABREVCOM.f
wget http://www.netlib.org/templates/double/STOPTEST2.f
wget http://www.netlib.org/lapack/util/dlamch.f

getbreak.fは拾い物． あとは，メインとなるプログラム，マトベク用プログラムのmatvec，前処理を解くpsolveを用意．

コンパイルはこんな感じで．blasが必要．

gfortran dlamch.f BiCGSTAB.f BiCGSTABREVCOM.f getbreak.f STOPTEST2.f main.f90 -lblas 

main.f90

module variables
  real(8),allocatable :: a(:,:),b(:),x(:),work(:,:)
  integer n
end module variables
 
module libmatrix
  use variables
contains
  subroutine MATVEC( ALPHA, X, BETA, Y )
    real(8),intent(in) :: alpha,beta
    real(8),intent(in) :: x(n)
    real(8),intent(out) :: y(n)
    y(:)=alpha*matmul(a(:,:),x(:))+beta*x(:)
  end subroutine MATVEC
 
  subroutine PSOLVE( X, B )
    real(8),intent(in) :: b(n)
    real(8),intent(out) :: x(n)
    x(:)=b(:)
 
  end subroutine PSOLVE
 
end module libmatrix
 
program main
  use libmatrix
  integer ldw,iter,info
  real(8) :: resid=1.d-6
  n=3
  ldw=3
  iter=3
 
  allocate(a(n,n),b(n),x(n),work(ldw,7))
  a(1,:)=(/1.d0,0.d0,0.d0/)
  a(2,:)=(/1.d0,1.d0,0.d0/)
  a(3,:)=(/1.d0,0.d0,1.d0/)
  b(:)=(/1.d0,2.d0,3.d0/)
  call BICGSTAB(3,B, X, WORK, LDW, ITER, RESID,MATVEC,PSOLVE, INFO )
  write(*,*) x
end program main

getbreak.f

      DOUBLE PRECISION FUNCTION GETBREAK()
*
*     Get breakdown parameter tolerance; for the test routine,
*     set to machine precision.
*
      DOUBLE PRECISION EPS, DLAMCH
*
      EPS = DLAMCH('EPS')
      GETBREAK = EPS**2
*
      RETURN
      END