====== BiCGSTAB ======
適当にソースコードを拾ってくる．
<code bash>
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
</code>
getbreak.fは拾い物．
あとは，メインとなるプログラム，マトベク用プログラムのmatvec，前処理を解くpsolveを用意．

コンパイルはこんな感じで．blasが必要．
<code bash>
gfortran dlamch.f BiCGSTAB.f BiCGSTABREVCOM.f getbreak.f STOPTEST2.f main.f90 -lblas 
</code>

main.f90
<code fortran>
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
</code>
getbreak.f
<code fortran>
      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
</code>