V3FIT
polysvd.f
1 !CALL polysvd(kpp,strial,ptrial,SIZE(strial),cfnt(1:kpp))
2  SUBROUTINE polysvd(nfitin,x,y,m,b)
3 ! USAGE SVD regression of input vectors x, y of length m
4 ! producess nfitin polynoimial coefficients b
5  USE stel_kinds
6  USE stel_constants
7  IMPLICIT NONE
8 C-----------------------------------------------
9 C D u m m y A r g u m e n t s
10 C-----------------------------------------------
11  INTEGER, INTENT(in) :: m, nfitin
12  REAL (rprec), DIMENSION(m), INTENT(in) :: x,y
13  REAL(rprec) ,DIMENSION(nfitin), INTENT(out) :: b
14 C-----------------------------------------------
15 C L o c a l V a r i a b l e s
16 C-----------------------------------------------
17  REAL (rprec) :: cutoff=1.e-7_dp
18  REAL (rprec), DIMENSION(:,:), ALLOCATABLE
19  1 :: amatrix, vv ,uu, wwd
20  REAL (rprec), DIMENSION(:) , ALLOCATABLE :: ww
21  REAL (rprec), DIMENSION(:), ALLOCATABLE, SAVE :: apar
22  INTEGER , DIMENSION(:) , ALLOCATABLE :: pwr
23  INTEGER :: i, j, n=11
24  n=nfitin
25  ALLOCATE(amatrix(m,n))
26  ALLOCATE(uu(m,n))
27  ALLOCATE(pwr(n))
28  ALLOCATE(apar(n))
29  ALLOCATE(wwd(n,n))
30  ALLOCATE(ww(n))
31  ALLOCATE(vv(n,n))
32  DO i=1,n
33  pwr(i)=i-1
34  ENDDO
35 ! make amatrix
36  DO i=1,m
37  DO j=1,n
38  IF(x(i).eq.0.) THEN
39  amatrix(i,j)=zero
40  ELSE
41  amatrix(i,j)=x(i)**pwr(j)
42  ENDIF
43  ENDDO
44  ENDDO
45  uu=amatrix
46  CALL svdcmp(uu,m,n,m,n,ww,vv) ! m rows > n columns
47  DO i=1,n
48  DO j=1,n
49  wwd(i,j)=0 ! reset after matmul
50  ENDDO
51  ENDDO
52  DO i=1,n
53  wwd(i,i)=1/ww(i)
54  ENDDO
55  DO i= 1,SIZE(ww)
56  IF(ww(i) .lt. cutoff) wwd(i,i)=0
57  ENDDO
58  apar=0.
59  apar=matmul(vv,matmul(wwd,matmul(transpose(uu),y)))
60  b(1:n)=apar(1:n)
61  DEALLOCATE(amatrix)
62  DEALLOCATE(uu)
63  DEALLOCATE(pwr)
64  DEALLOCATE(apar)
65  DEALLOCATE(wwd)
66  DEALLOCATE(ww)
67  DEALLOCATE(vv)
68  END SUBROUTINE polysvd
69 
70 ! PROGRAM driver
71 ! TESTED against IDL with this PARAMETER set. Results agree.
72 ! USE stel_kinds
73 ! USE stel_constants
74 ! INTEGER, PARAMETER :: ns=50, kpp=5
75 ! REAL(rprec), DIMENSION(ns) :: s, y
76 ! REAL(rprec), DIMENSION(kpp) :: b
77 ! DO i=1,ns
78 ! s(i)=one*(i-1)/(ns-1)
79 ! ENDDO
80 ! y=COS(s*pio2)**2
81 ! CALL polysvd(kpp,s,y,SIZE(s),b)
82 ! WRITE(6,109)b
83 ! WRITE(6,109)s
84 ! WRITE (6,109)y
85 !109 FORMAT(/,' =[',3(x,1pe14.6,1h,),'$')
86 ! END PROGRAM
87 
88