Télécharger la plaquette Cast3M


fr en

Télécharger xsrfc.eso

Retour à la liste

Numérotation des lignes : 

xsrfc

  1. C XSRFC     SOURCE    PV090527  23/01/27    21:16:11     11574          
  2.       subroutine xsrfc(x,y,z,nobs,coeff,errmax,errmoy)
  3. c source : https://people.math.sc.edu/Burkardt/f_src/qrsolv/qrsolv.html
  4.  
  5.         implicit real*8 (a-h,o-z)
  6.         implicit integer (i-n)
  7.  
  8.       integer ncoeff,n
  9.       parameter (ncoeff = 10,n=10)
  10.       integer nobs,i,ind,itask,jpvt(n),kr
  11.  
  12.  
  13. c  nobs : nombre de lignes de la table des points de base x,y,z
  14. c  ncoeff : nombre de coefficients à déterminer.
  15. c  resi : residue de l'approximation
  16.       real*8 errmax,resi(nobs)
  17.       real*8 zc(nobs),zmin,zmax
  18.       real*8 x(nobs),y(nobs),z(nobs),a(nobs,ncoeff),coeff(ncoeff)
  19.       real*8 z0(nobs),a0(nobs,ncoeff),errmoy,qraux(n),tol,work(n)
  20.  
  21. c sauvegarde des z initiaux 
  22.       z0=z  
  23.  
  24. c construct design matrix A
  25.       a(:,1) = 1.0d0 
  26.       a(:,2) = x
  27.       a(:,3) = y
  28.       a(:,4) = x**2
  29.       a(:,5) = y**2
  30.       a(:,6) = x*y
  31.       a(:,7) = x**2*y
  32.       a(:,8) = x*y**2
  33.       a(:,9) = x**3
  34.       a(:,10) = y**3
  35.  
  36.       a0=a
  37. c print*,a
  38. c  call fit(a,z,coeff) ! need a subroutine that solves A*coeff = Z in the least-squares sense
  39. c 
  40. c  write ( *, '(a)' ) ' '
  41. c   write ( *, '(a)' ) 'DQRLS_TEST'
  42. c   write ( *, '(a)' ) '  DQRLS solves a linear system A*coeff = z '
  43. c   write ( *, '(a)' ) '  in the least squares sense.'
  44.  
  45.       tol = 1.0D-8
  46.  
  47. c    call r8mat_print ( nobs, n, a, '  Coefficient matrix A:' )
  48. c  
  49. c    call r8vec_print ( nobs, z, '  Right hand side z:' )
  50. c  
  51. c    Solve least-squares problem.
  52.  
  53.       itask = 1
  54.       call dqrls (a,nobs,nobs,n,tol,kr,z,coeff,resi,jpvt,qraux,work,
  55.      #itask, ind )
  56.  
  57. c  Print results.
  58. c
  59. c  write ( *, '(a)' ) ' '
  60. c  write ( *, '(a,i4)' ) '  Error code =', ind
  61. c  write ( *, '(a,i4)' ) '  Estimated matrix rank =', kr
  62. c
  63. c  call r8vec_print ( n, coeff, '  Least squares solution coeff:' )
  64. c 
  65. c x0=coeff(1)
  66. c x1=coeff(2)
  67. c y1=coeff(3)
  68. c x2=coeff(4)
  69. c y2=coeff(5)
  70. c xy1=coeff(6)
  71. c x2y=coeff(7)
  72. c xy2=coeff(8)
  73. c x3=coeff(9)
  74. c y3=coeff(10)
  75.  
  76.  
  77. c  print*, "k0m(x,y)=",coeff(1),"+(",coeff(2),")*x+(",coeff(3),")*y+(",coeff(4),")*x**2+(",coeff(5),")*y**2+(",coeff(6),")*x*y"&
  78. c  ,"+(",coeff(7),")*x**2*y+(",coeff(8),")*x*y**2+(",coeff(9),")*x**3+(",coeff(10),")*y**3"
  79.  
  80.  
  81. c  call r8vec_print ( nobs, z, '  Residuals A*coeff-z' )
  82.  
  83.       zc=matmul(a0,coeff)
  84.       zmin=minval(z0)
  85.       zmax=maxval(z0)
  86.  
  87.       do i=1,nobs
  88.         if(zc(i).lt.zmin) then
  89.         zc(i)=zmin
  90.         else if(zc(i).gt.zmax) then
  91.         zc(i)=zmax
  92.         end if
  93.       end do
  94.  
  95. c print*, zc 
  96. c print*, z0 
  97. c read*
  98.  
  99.       errmax =sqrt(maxval((zc-z0)**2/z0**2))
  100.       errmoy= sqrt(sum((zc-z0)**2/z0**2))/nobs
  101.  
  102.  
  103. c erreur maximale
  104. c   errmax=maxval(resi/z)
  105. c   print*, errmax
  106.  
  107.       return
  108.  
  109.  
  110.  
  111.       end 
  112.  
  113.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales