![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to sqrsl.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / view2d / num|
Return to sqrsl.f CVS log | Up to [Research Unix] / researchv10no / cmd / view2d / num |
1.1 root 1: subroutine sqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
2: integer ldx,n,k,job,info
3: real x(ldx,1),qraux(1),y(1),qy(1),qty(1),b(1),rsd(1),xb(1)
4: c
5: c sqrsl applies the output of sqrdc to compute coordinate
6: c transformations, projections, and least squares solutions.
7: c for k .le. min(n,p), let xk be the matrix
8: c
9: c xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k)))
10: c
11: c formed from columnns jpvt(1), ... ,jpvt(k) of the original
12: c n x p matrix x that was input to sqrdc (if no pivoting was
13: c done, xk consists of the first k columns of x in their
14: c original order). sqrdc produces a factored orthogonal matrix q
15: c and an upper triangular matrix r such that
16: c
17: c xk = q * (r)
18: c (0)
19: c
20: c this information is contained in coded form in the arrays
21: c x and qraux.
22: c
23: c on entry
24: c
25: c x real(ldx,p).
26: c x contains the output of sqrdc.
27: c
28: c ldx integer.
29: c ldx is the leading dimension of the array x.
30: c
31: c n integer.
32: c n is the number of rows of the matrix xk. it must
33: c have the same value as n in sqrdc.
34: c
35: c k integer.
36: c k is the number of columns of the matrix xk. k
37: c must nnot be greater than min(n,p), where p is the
38: c same as in the calling sequence to sqrdc.
39: c
40: c qraux real(p).
41: c qraux contains the auxiliary output from sqrdc.
42: c
43: c y real(n)
44: c y contains an n-vector that is to be manipulated
45: c by sqrsl.
46: c
47: c job integer.
48: c job specifies what is to be computed. job has
49: c the decimal expansion abcde, with the following
50: c meaning.
51: c
52: c if a.ne.0, compute qy.
53: c if b,c,d, or e .ne. 0, compute qty.
54: c if c.ne.0, compute b.
55: c if d.ne.0, compute rsd.
56: c if e.ne.0, compute xb.
57: c
58: c note that a request to compute b, rsd, or xb
59: c automatically triggers the computation of qty, for
60: c which an array must be provided in the calling
61: c sequence.
62: c
63: c on return
64: c
65: c qy real(n).
66: c qy conntains q*y, if its computation has been
67: c requested.
68: c
69: c qty real(n).
70: c qty contains trans(q)*y, if its computation has
71: c been requested. here trans(q) is the
72: c transpose of the matrix q.
73: c
74: c b real(k)
75: c b contains the solution of the least squares problem
76: c
77: c minimize norm2(y - xk*b),
78: c
79: c if its computation has been requested. (note that
80: c if pivoting was requested in sqrdc, the j-th
81: c component of b will be associated with column jpvt(j)
82: c of the original matrix x that was input into sqrdc.)
83: c
84: c rsd real(n).
85: c rsd contains the least squares residual y - xk*b,
86: c if its computation has been requested. rsd is
87: c also the orthogonal projection of y onto the
88: c orthogonal complement of the column space of xk.
89: c
90: c xb real(n).
91: c xb contains the least squares approximation xk*b,
92: c if its computation has been requested. xb is also
93: c the orthogonal projection of y onto the column space
94: c of x.
95: c
96: c info integer.
97: c info is zero unless the computation of b has
98: c been requested and r is exactly singular. in
99: c this case, info is the index of the first zero
100: c diagonal element of r and b is left unaltered.
101: c
102: c the parameters qy, qty, b, rsd, and xb are not referenced
103: c if their computation is not requested and in this case
104: c can be replaced by dummy variables in the calling program.
105: c to save storage, the user may in some cases use the same
106: c array for different parameters in the calling sequence. a
107: c frequently occuring example is when one wishes to compute
108: c any of b, rsd, or xb and does not need y or qty. in this
109: c case one may identify y, qty, and one of b, rsd, or xb, while
110: c providing separate arrays for anything else that is to be
111: c computed. thus the calling sequence
112: c
113: c call sqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info)
114: c
115: c will result in the computation of b and rsd, with rsd
116: c overwriting y. more generally, each item in the following
117: c list contains groups of permissible identifications for
118: c a single callinng sequence.
119: c
120: c 1. (y,qty,b) (rsd) (xb) (qy)
121: c
122: c 2. (y,qty,rsd) (b) (xb) (qy)
123: c
124: c 3. (y,qty,xb) (b) (rsd) (qy)
125: c
126: c 4. (y,qy) (qty,b) (rsd) (xb)
127: c
128: c 5. (y,qy) (qty,rsd) (b) (xb)
129: c
130: c 6. (y,qy) (qty,xb) (b) (rsd)
131: c
132: c in any group the value returned in the array allocated to
133: c the group corresponds to the last member of the group.
134: c
135: c linpack. this version dated 08/14/78 .
136: c g.w. stewart, university of maryland, argonne national lab.
137: c
138: c sqrsl uses the following functions and subprograms.
139: c
140: c blas saxpy,scopy,sdot
141: c fortran abs,min0,mod
142: c
143: c internal variables
144: c
145: integer i,j,jj,ju,kp1
146: real sdot,t,temp
147: logical cb,cqy,cqty,cr,cxb
148: c
149: c
150: c set info flag.
151: c
152: info = 0
153: c
154: c determine what is to be computed.
155: c
156: cqy = job/10000 .ne. 0
157: cqty = mod(job,10000) .ne. 0
158: cb = mod(job,1000)/100 .ne. 0
159: cr = mod(job,100)/10 .ne. 0
160: cxb = mod(job,10) .ne. 0
161: ju = min0(k,n-1)
162: c
163: c special action when n=1.
164: c
165: if (ju .ne. 0) go to 40
166: if (cqy) qy(1) = y(1)
167: if (cqty) qty(1) = y(1)
168: if (cxb) xb(1) = y(1)
169: if (.not.cb) go to 30
170: if (x(1,1) .ne. 0.0e0) go to 10
171: info = 1
172: go to 20
173: 10 continue
174: b(1) = y(1)/x(1,1)
175: 20 continue
176: 30 continue
177: if (cr) rsd(1) = 0.0e0
178: go to 250
179: 40 continue
180: c
181: c set up to compute qy or qty.
182: c
183: if (cqy) call scopy(n,y,1,qy,1)
184: if (cqty) call scopy(n,y,1,qty,1)
185: if (.not.cqy) go to 70
186: c
187: c compute qy.
188: c
189: do 60 jj = 1, ju
190: j = ju - jj + 1
191: if (qraux(j) .eq. 0.0e0) go to 50
192: temp = x(j,j)
193: x(j,j) = qraux(j)
194: t = -sdot(n-j+1,x(j,j),1,qy(j),1)/x(j,j)
195: call saxpy(n-j+1,t,x(j,j),1,qy(j),1)
196: x(j,j) = temp
197: 50 continue
198: 60 continue
199: 70 continue
200: if (.not.cqty) go to 100
201: c
202: c compute trans(q)*y.
203: c
204: do 90 j = 1, ju
205: if (qraux(j) .eq. 0.0e0) go to 80
206: temp = x(j,j)
207: x(j,j) = qraux(j)
208: t = -sdot(n-j+1,x(j,j),1,qty(j),1)/x(j,j)
209: call saxpy(n-j+1,t,x(j,j),1,qty(j),1)
210: x(j,j) = temp
211: 80 continue
212: 90 continue
213: 100 continue
214: c
215: c set up to compute b, rsd, or xb.
216: c
217: if (cb) call scopy(k,qty,1,b,1)
218: kp1 = k + 1
219: if (cxb) call scopy(k,qty,1,xb,1)
220: if (cr .and. k .lt. n) call scopy(n-k,qty(kp1),1,rsd(kp1),1)
221: if (.not.cxb .or. kp1 .gt. n) go to 120
222: do 110 i = kp1, n
223: xb(i) = 0.0e0
224: 110 continue
225: 120 continue
226: if (.not.cr) go to 140
227: do 130 i = 1, k
228: rsd(i) = 0.0e0
229: 130 continue
230: 140 continue
231: if (.not.cb) go to 190
232: c
233: c compute b.
234: c
235: do 170 jj = 1, k
236: j = k - jj + 1
237: if (x(j,j) .ne. 0.0e0) go to 150
238: info = j
239: c ......exit
240: go to 180
241: 150 continue
242: b(j) = b(j)/x(j,j)
243: if (j .eq. 1) go to 160
244: t = -b(j)
245: call saxpy(j-1,t,x(1,j),1,b,1)
246: 160 continue
247: 170 continue
248: 180 continue
249: 190 continue
250: if (.not.cr .and. .not.cxb) go to 240
251: c
252: c compute rsd or xb as required.
253: c
254: do 230 jj = 1, ju
255: j = ju - jj + 1
256: if (qraux(j) .eq. 0.0e0) go to 220
257: temp = x(j,j)
258: x(j,j) = qraux(j)
259: if (.not.cr) go to 200
260: t = -sdot(n-j+1,x(j,j),1,rsd(j),1)/x(j,j)
261: call saxpy(n-j+1,t,x(j,j),1,rsd(j),1)
262: 200 continue
263: if (.not.cxb) go to 210
264: t = -sdot(n-j+1,x(j,j),1,xb(j),1)/x(j,j)
265: call saxpy(n-j+1,t,x(j,j),1,xb(j),1)
266: 210 continue
267: x(j,j) = temp
268: 220 continue
269: 230 continue
270: 240 continue
271: 250 continue
272: return
273: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.