| changeset 17 | 4b7486d77b39 |
| parent 0 | a5233793bf69 |
| 15:ebd8a4ddd7f2 | 17:4b7486d77b39 |
|---|---|
1 #====================================================================== |
1 #====================================================================== |
2 # S V D C M P . P L |
2 # S V D C M P . P L |
3 # doc: Sun Aug 1 09:51:37 1999 |
3 # doc: Sun Aug 1 09:51:37 1999 |
4 # dlm: Thu Jul 19 09:45:52 2001 |
4 # dlm: Wed Mar 11 16:54:26 2015 |
5 # (c) 1999 A.M. Thurnherr |
5 # (c) 1999 A.M. Thurnherr |
6 # uE-Info: 184 18 NIL 0 0 72 2 2 4 NIL ofnI |
6 # uE-Info: 188 1 NIL 0 0 72 10 2 4 NIL ofnI |
7 #====================================================================== |
7 #====================================================================== |
8 |
8 |
9 # SVDCMP routine from Numerical Recipes adapted to ANTS |
9 # SVDCMP routine from Numerical Recipes adapted to ANTS |
10 |
10 |
11 # HISTORY: |
11 # HISTORY: |
12 # Aug 01, 1999: - manually converted from c-source |
12 # Aug 01, 1999: - manually converted from c-source |
13 # Jul 19, 2001: - done *something* (message only?) |
|
14 # Mar 11, 2015: - picked up this never-used code because of need to fit circles |
|
15 |
|
13 |
16 |
14 # Notes: |
17 # Notes: |
15 # - everything passed as refs |
18 # - everything passed as refs |
16 |
19 |
17 require "$ANTS/nrutil.pl"; |
20 require "$ANTS/nrutil.pl"; |
18 require "$ANTS/pythag.pl"; |
21 require "$ANTS/pythag.pl"; |
19 |
22 |
23 use strict; |
|
24 |
|
20 sub svdcmp($$$) |
25 sub svdcmp($$$) |
21 { |
26 { |
22 my($aR,$wR,$vR) = @_; # params |
27 my($aR,$wR,$vR) = @_; # params |
23 my($flag,$i,$its,$j,$jj,$k,$l,$nm); # int |
28 my($flag,$i,$its,$j,$jj,$k,$l,$nm); # int |
24 my($anorm,$c,$f,$g,$h,$s,$scale,$x,$y,$z); # float |
29 my($anorm,$c,$f,$g,$h,$s,$scale,$x,$y,$z); # float |
25 my(@rv1); # float[] |
30 my(@rv1); # float[] |
26 |
31 |
27 vector(\@rv1,1,$#{$vR}); |
32 vector(\@rv1,1,$#{$vR}); |
33 $g = $scale = $anorm = 0; |
|
28 for ($i=1; $i<=$#{$vR}; $i++) { |
34 for ($i=1; $i<=$#{$vR}; $i++) { |
29 $l = $i+1; |
35 $l = $i+1; |
30 $rv1[$i] = $scale*$g; |
36 $rv1[$i] = $scale*$g; |
31 $g = 0; $s = 0; $scale = 0; |
37 $g = $s = $scale = 0; |
32 if ($i <= $#{$aR}) { |
38 if ($i <= $#{$aR}) { |
33 for ($k=$i; $k<=$#{$aR}; $k++) { |
39 for ($k=$i; $k<=$#{$aR}; $k++) { |
34 $scale += abs($aR->[$k][$i]); |
40 $scale += abs($aR->[$k][$i]); |
35 } |
41 } |
36 if ($scale) { |
42 if ($scale) { |
85 for ($k=$l; $k<=$#{$vR}; $k++) { |
91 for ($k=$l; $k<=$#{$vR}; $k++) { |
86 $aR->[$i][$k] *= $scale; |
92 $aR->[$i][$k] *= $scale; |
87 } |
93 } |
88 } |
94 } |
89 } |
95 } |
90 $anorm = &FMAX($anorm,(abs($wR->[$i])+abs($rv1[$i]))); |
96 $anorm = MAX($anorm,(abs($wR->[$i])+abs($rv1[$i]))); |
91 } |
97 } |
92 for ($i=$#{$vR}; $i>=1; $i--) { |
98 for ($i=$#{$vR}; $i>=1; $i--) { |
93 if ($i < $#{$vR}) { |
99 if ($i < $#{$vR}) { |
94 if ($g) { |
100 if ($g) { |
95 for ($j=$l; $j<=$#{$vR}; $j++) { |
101 for ($j=$l; $j<=$#{$vR}; $j++) { |
102 for ($k=$l; $k<=$#{$vR}; $k++) { |
108 for ($k=$l; $k<=$#{$vR}; $k++) { |
103 $vR->[$k][$j] += $s*$vR->[$k][$i]; |
109 $vR->[$k][$j] += $s*$vR->[$k][$i]; |
104 } |
110 } |
105 } |
111 } |
106 } |
112 } |
107 for ($j=$l; $j<=$#{$vR; $j++) { |
113 for ($j=$l; $j<=$#{$vR}; $j++) { |
108 $vR->[$i][$j] = 0; $vR->[$j][$i] = 0; |
114 $vR->[$i][$j] = 0; $vR->[$j][$i] = 0; |
109 } |
115 } |
110 } |
116 } |
111 $vR->[$i][$i] = 1; |
117 $vR->[$i][$i] = 1; |
112 $g = $rv1[$i]; |
118 $g = $rv1[$i]; |
113 $l = $i; |
119 $l = $i; |
114 } |
120 } |
115 for ($i=IMIN($#{$aR},$#{$vR}); $i>=1; $i--) { |
121 for ($i=MIN($#{$aR},$#{$vR}); $i>=1; $i--) { |
116 $l = $i+1; |
122 $l = $i+1; |
117 $g = $wR->[$i]; |
123 $g = $wR->[$i]; |
118 for ($j=$l; $j<=$#{$vR}; $j++) { |
124 for ($j=$l; $j<=$#{$vR}; $j++) { |
119 $aR->[$i][$j] = 0; |
125 $aR->[$i][$j] = 0; |
120 } |
126 } |
144 $flag = 1; |
150 $flag = 1; |
145 for ($l=$k; $l>=1; $l--) { |
151 for ($l=$k; $l>=1; $l--) { |
146 $nm = $l-1; |
152 $nm = $l-1; |
147 if ((abs($rv1[$l])+$anorm) == $anorm) { |
153 if ((abs($rv1[$l])+$anorm) == $anorm) { |
148 $flag = 0; |
154 $flag = 0; |
149 break; |
155 last; |
150 } |
156 } |
151 break if ((abs($wR->[$nm])+$anorm) == $anorm); |
157 last if ($nm == 0) || ((abs($wR->[$nm])+$anorm) == $anorm); ## $nm == 0 test not in original code |
152 } |
158 } |
153 if ($flag) { |
159 if ($flag) { |
154 $c = 0; |
160 $c = 0; |
155 $s = 1; |
161 $s = 1; |
156 for ($i=$l; $i<=$k; $i++) { |
162 for ($i=$l; $i<=$k; $i++) { |
157 $f = $s*$rv1[$i]; |
163 $f = $s*$rv1[$i]; |
158 $rv1[$i] = $c*$rv1[$i]; |
164 $rv1[$i] = $c*$rv1[$i]; |
159 break if ((abs($f)+$anorm) == $anorm); |
165 last if ((abs($f)+$anorm) == $anorm); |
160 $g = $wR->[$i]; |
166 $g = $wR->[$i]; |
161 $h = &pythag($f,$g); |
167 $h = &pythag($f,$g); |
162 $wR->[$i] = $h; |
168 $wR->[$i] = $h; |
163 $h = 1/$h; |
169 $h = 1/$h; |
164 $c = $g*$h; |
170 $c = $g*$h; |
177 $wR->[$k] = -$z; |
183 $wR->[$k] = -$z; |
178 for ($j=1; $j<=$#{$vR}; $j++) { |
184 for ($j=1; $j<=$#{$vR}; $j++) { |
179 $vR->[$j][$k] = -$vR->[$j][$k]; |
185 $vR->[$j][$k] = -$vR->[$j][$k]; |
180 } |
186 } |
181 } |
187 } |
182 break; |
188 # print(STDERR "its($k) = $its\n"); |
189 last; |
|
183 } |
190 } |
184 croak("no convergence in 30 svdcmp iterations\n") if ($its == 30); |
191 croak("no convergence in 30 svdcmp iterations\n") if ($its == 30); |
185 $x = $wR->[$l]; |
192 $x = $wR->[$l]; |
186 $nm = $k-1; |
193 $nm = $k-1; |
187 $y = $wR->[$nm]; |
194 $y = $wR->[$nm]; |