用FORTRAN怎样对矩阵做LU分解

huangfen1234

请问各位，用FORTRAN怎样对矩阵做LU分解，我在IMSL里面好像没看到

swordlj

不要局限于某种语言，算法是关键！

huangfen1234

不是的，我现在是F程序里面要用到，想直接从IMSL调用这个功能的函数，在一本算法集上面看到应该是LUDCMP函数，但IMSL里搜不到，不知道是不是盗版软件的原因，难道是IMSL不全，要自己把子程序编了也可以，但没那个必要

lishhit

imsl里面有 不过请尽量支持正版

lishhit

在何光渝的《VISUAL FORTRAN 常用数值算法集》的第一章里也有，现成的一个源程序。

huangfen1234

呵呵，谢谢各位了，我的矩阵是复数矩阵，所以麻烦一点，现在LU已经解决，还要用到QR分解，IMSL里的LQRRR只对实矩阵分解，不像LU有另外的函数对复矩阵分解，问题还没解决，有哪位知道的话指点一下，感激不尽，谢谢。

zxb0412

到网上下载一个 Blas 和 Lapack 库（www.netlib.org），有针对不同 CPU 及相应平台下二进制文件，当然也有 Windows 下的函数库。
个人觉得 Blas 库或 Lapack 库的稳定性和速度远优于算法集中的函数。你所提到的几种分解，在 Lapack 中都有。

andytsai

huangfen,
你可以看看IMSL中的ORTH函数，是否可以解决你的问题。
或是直接到IMSL的论坛
http://www.vni.com.tw/forums/
发问~

[ 本帖最后由 andytsai 于 2008-4-29 12:27 编辑 ]

piteqiu

长见识了，好多库都不懂

ma

高斯分解
SUBROUTINE sparin_gauss(kv,kdiag)
!
! This subroutine performs Gaussian factorisation of a skyline matrix.
!
IMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER,INTENT(IN)::kdiag(
REAL(iwp),INTENT(OUT)::kv(:)
REAL(iwp)::num,den,fac,zero=0.0_iwp
INTEGER::n,ii,i,j,k,l,kk,l1,l2,l3
n=UBOUND(kdiag,1)
DO j=1,n-1
   den=kv(kdiag(j))
   ii=0                 
   DO i=j+1,n
     ii=ii+1
     l=kdiag(i)-ii
     IF(l-kdiag(i-1)>zero)THEN
       num=kv(l)
       fac=num/den
       kk=-1
       DO k=i,n
         kk=kk+1
         l1=kdiag(i+kk)-kk
         l2=l1-ii
         l3=kdiag(i+kk-1)
         IF(l2-l3>zero)kv(l1)=kv(l1)-fac*kv(l2)
       END DO
     END IF
   END DO
END DO
RETURN
END SUBROUTINE sparin_gauss

ma

subroutine banred(bk,n)
! gaussian reduction on a vector stored as an upper triangle
implicit none
real,intent(in out)::bk(;integer,intent(in)::n
integer::i,il1,kbl,j,ij,nkb,m,ni,nj,iw ; real::sum
iw = ubound(bk,1)/n-1
       do i=2,n
          il1=i-1;kbl=il1+iw+1
          if(kbl-n>0)kbl=n
          do j=i,kbl
             ij=(j-i)*n+i;sum=bk(ij);nkb=j-iw
             if(nkb<=0)nkb=1
             if(nkb-il1<=0)then
                do m=nkb,il1
                   ni=(i-m)*n+m ; nj=(j-m)*n+m
                   sum=sum-bk(ni)*bk(nj)/bk(m)
                end do
             end if
             bk(ij)=sum
           end do
       end do
return
end subroutine banred

fhuang

谢谢各位！