Fortran Coder

标题: 试试 高亮 [打印本页]

作者: 珊瑚虫    时间: 2014-1-23 22:05
标题: 试试 高亮
[Fortran] 纯文本查看 复制代码
module lsqcurvefit
!----------------------------------------module coment
!  Version     :  V1.0   
!  Coded by    :  syz
!  Date        :  2010.07.09
!-----------------------------------------------------
!  Description :   任意阶多项式曲线拟合模块
!  
!  Post Script :
!      1.
!      2.
!
!-----------------------------------------------------
!  Contains    :
!      1.    拟合函数,基函数
!      2.    最小二乘法函数
!-----------------------------------------------------
!  Parameters  :
!      1.
!      2.
!-----------------------------------------------------
contains
subroutine solve(x,y,N,c,m)
!---------------------------------subroutine  comment
!  Version   :  V1.0   
!  Coded by  :  syz
!  Date      :  
!-----------------------------------------------------
!  Purpose     :  任意阶多项式拟合函数
!
!  Post Script :
!       1.
!       2.
!       3.   
!-----------------------------------------------------
!  Input  parameters  :
!       1.  x,y --- 输入数据
!       2.  N --- x,y向量的维数
!       3.  m----希望用m阶多项式拟合
!  Output parameters  :
!       1.
!       2.
!  Common parameters  :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)

integer::N,m
real*8::x(n),y(n),c(M)

real*8::bv(M)

integer::i

!A系数矩阵
real*8::A(N,M)


!如果多项式阶数加1高于实测数据个数则报错
!-------------
if (m>n) then

write(11,101)
stop
end if

101 format(/,'warning:the order of polynomial+1 > the ',/,&
              'number of date,that is forbidden')
!------------------------------



do i=1,n
   
    call basefunc(bv,x(i),m)
   
    A(i,:)=bv

end do

call leas_eq(A,y,c,N,m)


end subroutine solve



subroutine basefunc(bv,x,m)
!---------------------------------subroutine  comment
!  Version   :  V1.0   
!  Coded by  :  syz
!  Date      :2010.07.09  
!-----------------------------------------------------
!  Purpose     :  任意阶多项式基函数
!
!  Post Script :
!       1.
!       2.
!       3.   
!-----------------------------------------------------
!  Input  parameters  :
!       1.    x 标量
!       2.    m 多项式阶数
!  Output parameters  :
!       1.
!       2.
!  Common parameters  :
!       1.
!       2.
!----------------------------------------------------
implicit real*8(a-z)

integer::M

real*8::bv(m)

integer::i

bv(1)=1d0

do  i=2,m
    bv(i)=bv(i-1)*x
end do
end subroutine basefunc

subroutine  leas_eq(A,b,x,M,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0   
!  Coded by  :  syz
!  Date      :  
!-----------------------------------------------------
!  Purpose   :  通过修正的Gram-Schmidt正交化求最小二问题
!               方法函数
!-----------------------------------------------------
!  Method    :
!               对超定方程 A进行QR分解后  方程变为
!                   QR x=b
!                => Rx=Q'b   R为上三角阵
!                => 回带,可以求得最小二乘意义下的解
!-----------------------------------------------------
!  Post Script :
!       1.       即求解超定方程组 Ax=b    其中A(M,N)  M>N   
!       2.
!----------------------------------------------------
implicit real*8(a-z)

integer::M,N
real*8::A(M,N),Q(M,N),R(N,N)
real*8::b(M)
real*8::QT(N,M)  !Q的转置矩阵
real*8::QTb(N)   !Q'b
real*8::x(N)

call gram_dec(A,Q,R,M,N)

QT=transpose(Q)
QTb=matmul(QT,b)  !  Rx=Q'b

call uptri(R,QTb,x,N) !回带

end subroutine leas_eq

subroutine gram_dec(A,Q,R,M,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0   
!  Coded by  :  syz
!  Date      :  
!-----------------------------------------------------
!  Purpose   :   采用修正的Gram-Schmidt分解求矩阵的QR分解
!   
!-----------------------------------------------------
!  Input  parameters  :
!       1.    A原始矩阵
!       2.    A(M,N)
!  Output parameters  :
!       1.    分解结果为  Q(M,N):注意Q不是方阵,Q列向量为标准正交基
!       2.                 R(N,N):R是方阵
!       3.   
!----------------------------------------------------
!  Post Script :
!       1.  注意矩阵的维数,分解后Q列向量是正交的
!       2.  关于编程方法可以参看《矩阵分析与应用》张贤达编著
!       3.  详细的数学解释,可以参看麻省理工学院的
!           线性代数教材《Linear Algebra with Application》
!----------------------------------------------------
implicit real*8(a-z)

integer::M,N
integer::i,j,k

real*8::A(M,N),Q(M,N),R(N,N)

real*8::vec_temp(M)



R(1,1)=dsqrt(dot_product(a(:,1),a(:,1)))

Q(:,1)=a(:,1)/R(1,1)


do k=2,N

      do j=1,k-1
        R(j,k)=dot_product(Q(:,j),A(:,k))   
      end do
   
      vec_temp=A(:,k)
   
      do j=1,k-1
   
        vec_temp=vec_temp-Q(:,j)*R(j,k)
   
      end do


    R(k,k)=dsqrt(dot_product(vec_temp,vec_temp))

     
    Q(:,k)=vec_temp/R(k,k)

end do

end subroutine gram_dec

subroutine uptri(A,b,x,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0   
!  Coded by  :  syz
!  Date      :  2010-4-8
!-----------------------------------------------------
!  Purpose   :  上三角方程组的回带方法
!                 Ax=b
!-----------------------------------------------------
!  Input  parameters  :
!       1.   A(N,N)系数矩阵
!       2.   b(N)右向量
!       3.   N方程维数
!  Output parameters  :
!       1.  x  方程的根
!       2.
!  Common parameters  :
!
!----------------------------------------------------

implicit real*8(a-z)

integer::i,j,k,N

real*8::A(N,N),b(N),x(N)

x(N)=b(N)/A(N,N)

!回带部分
do i=n-1,1,-1
   
    x(i)=b(i)
   do j=i+1,N
    x(i)=x(i)-a(i,j)*x(j)
   end do
    x(i)=x(i)/A(i,i)

end do

end subroutine uptri

end module lsqcurvefit



program main
!--------------------------------------program comment
!  Version   :  V1.0   
!  Coded by  :  syz
!  Date      :  2010.07.09
!-----------------------------------------------------
!  Purpose     : 多项式拟合主函数
!
!  Post Script :
!       1.
!       2.
!   
!-----------------------------------------------------
!  In put data  files :
!       1.
!       2.
!  Output data files  :
!       1.
!       2.
!
!-----------------------------------------------------
use  lsqcurvefit

implicit real*8(a-z)

real*8::x(7),y(7)

real*8::c1(3),c2(4),c3(9)

open(unit=11,file='result.txt')

x=(/-3d0,-2d0,-1d0,0d0,1d0,2d0,3d0/)
y=(/4d0,2d0,3d0,0d0,-1d0,-2d0,-5d0/)

call solve(x,y,7,c1,3)

call solve(x,y,7,c2,4)



write(11,101)c1,c2

101 format(/T4,'多项式拟合曲线拟合',//,&
              T3,'采用二阶多项式拟合,系数为:',3(/,F16.11),//,&
              T3,'采用三阶多项式拟合,系数为:',4(/,F16.11),/)
              

call solve(x,y,7,c3,9)              
              

end program main



作者: fcode    时间: 2014-1-23 22:08
看起来效果还可以哦
作者: 珊瑚虫    时间: 2014-1-23 22:08
8错8错 很不错
作者: 珊瑚虫    时间: 2014-1-23 22:14
fcode 发表于 2014-1-23 22:08
看起来效果还可以哦

固定格式呢?
作者: fcode    时间: 2014-1-23 22:18
珊瑚虫 发表于 2014-1-23 22:14
固定格式呢?

你可以试试,不敢说100%,大多数还是支持
作者: Colderheart    时间: 2014-1-23 22:19
我来顶顶,灌灌水:P:P:P:P:P
作者: 珊瑚虫    时间: 2014-1-23 22:24
[Fortran] 纯文本查看 复制代码
      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SX(1),SY(1),STEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
C
      STEMP = 0.0E0
      SDOT = 0.0E0
      IF(N.LE.0)RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF(INCX.LT.0)IX = (-N+1)*INCX + 1
      IF(INCY.LT.0)IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        STEMP = STEMP + SX(IX)*SY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      SDOT = STEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF( M .EQ. 0 ) GO TO 40
      DO 30 I = 1,M
        STEMP = STEMP + SX(I)*SY(I)
   30 CONTINUE
      IF( N .LT. 5 ) GO TO 60
   40 MP1 = M + 1
      DO 50 I = MP1,N,5
        STEMP = STEMP + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) +
     *   SX(I + 2)*SY(I + 2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4)
   50 CONTINUE
   60 SDOT = STEMP
      RETURN
      END

作者: jeydragon    时间: 2014-5-12 22:40
[Fortran] 纯文本查看 复制代码
program main
implicit none
end program main


作者: jeydragon    时间: 2014-5-12 22:41
总算成功了,原来是在高级模式里面。
作者: 楚香饭    时间: 2014-5-12 22:43
jeydragon 发表于 2014-5-12 22:41
总算成功了,原来是在高级模式里面。

<> 这个工具就可以,普通模式也有的哦
作者: jeydragon    时间: 2014-5-12 23:09
<program example>
<implicit none>
<write(*,*) 'Hello, World'>
<end program example>
作者: fcode    时间: 2014-5-12 23:17
你没看到这些东西?其中有一个 <> 图标

作者: jason388    时间: 2014-11-5 08:20
本帖最后由 jason388 于 2014-11-5 08:21 编辑

[Fortran] 纯文本查看 复制代码
print*,'Hello Fortran!';end





欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2