! How to use
!  real*8 :: c(n_dim,n_dim),e(n_dim,n_dim)
!  integer :: n=4
!     n           size of the matrix to be inversed
!     n_dim       array size
!     e           unit matrix
!     c           matrix to be inversed
!     call gaussr_inverse(c,n,n_dim,e,n,n_dim)

!flag
subroutine gaussr_inverse(a,n,np,b,m,mp)
  !     Commented, reorganized and understood by xxxx, 1999-11-03
  !     1.   very good designed scheme.
  !          seems No room for performance improve
  !     2.   row normalization (implicite pivoting) has not much  effect
  !           for full pivoting procedure

  !     991112  ASY/FJT, xxxx
  !     Test again. even important for full pivoting procedure

  !     3.   row normalization (implicite pivoting) has significant
  !           accuracy enhancement on partial pivoting procedure
  !              7.e-5 -> 2.e-11

  !     000623  ASY/FJT, xxx
  !     Real version


  implicit none
  integer :: n,np,m,mp,i,j,k,irow,icol,l
  real*8 a(np,np),b(np,mp),pivinv,cc
  real*8, allocatable, dimension(:) :: dum
  real*8 big
  integer, allocatable, dimension(:) :: ipiv,indxr,indxc
  real*8, allocatable, dimension(:) :: c1

  allocate(ipiv(1:n),indxr(1:n),indxc(1:n),dum(1:max(n,m)),c1(1:n))
  ipiv=0

  !     row normalization (implicite pivoting)
  do i=1,n
     c1(i)=1/sqrt(sum(abs(a(i,1:n))**2))
     a(i,1:n)=a(i,1:n)*c1(i)
     b(i,1:m)=b(i,1:m)*c1(i)
  end do

  do i=1,n

     !     Full pivoting, however, not necessarily starting from 1st col
     big=0.d0
     do j=1,n
        if(ipiv(j)==1)cycle
        do k=1,n
           if (ipiv(k)==1) cycle
           if (abs(a(j,k)) <= big) cycle
           big=abs(a(j,k))
           irow=j
           icol=k
        end do
     end do

     !     ipiv(icol)=1 means the column icol has been selected once
     ipiv(icol)=ipiv(icol)+1
     if (irow.ne.icol) then

        !     exchange a(irow,:) with a(icol,:)
        !     so a(irow,irow) is the selected pivot
        dum(1:n)= a(irow,1:n)  
        a(irow,1:n)= a(icol,1:n)  
        a(icol,1:n)=dum(1:n)

        dum(1:m)= b(irow,1:m)  
        b(irow,1:m)= b(icol,1:m)  
        b(icol,1:m)=dum(1:m)
     endif

     !     bookkeeping row- and column-indices
     indxr(i)=irow
     indxc(i)=icol
     if (abs(a(icol,icol))==0.d0) pause 'singular matrix.'
     pivinv=1.d0/a(icol,icol)
     a(icol,icol)=1.d0
     a(icol,1:n)=a(icol,1:n)*pivinv
     b(icol,1:m)=b(icol,1:m)*pivinv

     do l=1,n
        if(l==icol)cycle  ! only for row l not equal to icol
        cc=a(l,icol)
        a(l,icol)=0.d0    ! clever design
        a(l,1:n)=a(l,1:n)-a(icol,1:n)*cc
        b(l,1:m)=b(l,1:m)-b(icol,1:m)*cc
     end do
  end do

  do l=1,n
     if(indxr(l) == indxc(l)) cycle
     dum(1:n)=a(1:n,indxr(l))
     a(1:n,indxr(l))=a(1:n,indxc(l))
     a(1:n,indxc(l))=dum(1:n)
  end do

  do i=1,n
     a(1:n,i)=a(1:n,i)*c1(i)
  end do

  deallocate(ipiv,indxr,indxc,dum,c1)
end subroutine gaussr_inverse

2楼: Originally posted by pippi6 at 2013-06-05 15:54:26
给你个程序

! How to use
!  real*8 :: c(n_dim,n_dim),e(n_dim,n_dim)
!  integer :: n=4
!     n           size of the matrix to be inversed
!     n_dim       array size
!     e           unit ma ...

4楼: Originally posted by pippi6 at 2013-06-06 17:00:51
你这两句是想干什么？
                       ImCmplx(Row,j)=(Im(j),0)
                         ImCmplx(Row,j)=(0,0)
如果是想表示一个复数由实部和虚部构成，应这么写
                      ImCmplx(Row ...

4楼: Originally posted by pippi6 at 2013-06-06 17:00:51
你这两句是想干什么？
                       ImCmplx(Row,j)=(Im(j),0)
                         ImCmplx(Row,j)=(0,0)
如果是想表示一个复数由实部和虚部构成，应这么写
                      ImCmplx(Row ...