![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to sqrdc.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 sqrdc.f CVS log | Up to [Research Unix] / researchv10no / cmd / view2d / num |
1.1 root 1: subroutine sqrdc(x,ldx,n,p,qraux,jpvt,work,job)
2: integer ldx,n,p,job
3: integer jpvt(1)
4: real x(ldx,1),qraux(1),work(1)
5: c
6: c sqrdc uses householder transformations to compute the qr
7: c factorization of an n by p matrix x. column pivoting
8: c based on the 2-norms of the reduced columns may be
9: c performed at the users option.
10: c
11: c on entry
12: c
13: c x real(ldx,p), where ldx .ge. n.
14: c x contains the matrix whose decomposition is to be
15: c computed.
16: c
17: c ldx integer.
18: c ldx is the leading dimension of the array x.
19: c
20: c n integer.
21: c n is the number of rows of the matrix x.
22: c
23: c p integer.
24: c p is the number of columns of the matrix x.
25: c
26: c jpvt integer(p).
27: c jpvt contains integers that control the selection
28: c of the pivot columns. the k-th column x(k) of x
29: c is placed in one of three classes according to the
30: c value of jpvt(k).
31: c
32: c if jpvt(k) .gt. 0, then x(k) is an initial
33: c column.
34: c
35: c if jpvt(k) .eq. 0, then x(k) is a free column.
36: c
37: c if jpvt(k) .lt. 0, then x(k) is a final column.
38: c
39: c before the decomposition is computed, initial columns
40: c are moved to the beginning of the array x and final
41: c columns to the end. both initial and final columns
42: c are frozen in place during the computation and only
43: c free columns are moved. at the k-th stage of the
44: c reduction, if x(k) is occupied by a free column
45: c it is interchanged with the free column of largest
46: c reduced norm. jpvt is not referenced if
47: c job .eq. 0.
48: c
49: c work real(p).
50: c work is a work array. work is not referenced if
51: c job .eq. 0.
52: c
53: c job integer.
54: c job is an integer that initiates column pivoting.
55: c if job .eq. 0, no pivoting is done.
56: c if job .ne. 0, pivoting is done.
57: c
58: c on return
59: c
60: c x x contains in its upper triangle the upper
61: c triangular matrix r of the qr factorization.
62: c below its diagonal x contains information from
63: c which the orthogonal part of the decomposition
64: c can be recovered. note that if pivoting has
65: c been requested, the decomposition is not that
66: c of the original matrix x but that of x
67: c with its columns permuted as described by jpvt.
68: c
69: c qraux real(p).
70: c qraux contains further information required to recover
71: c the orthogonal part of the decomposition.
72: c
73: c jpvt jpvt(k) contains the index of the column of the
74: c original matrix that has been interchanged into
75: c the k-th column, if pivoting was requested.
76: c
77: c linpack. this version dated 08/14/78 .
78: c g.w. stewart, university of maryland, argonne national lab.
79: c
80: c sqrdc uses the following functions and subprograms.
81: c
82: c blas saxpy,sdot,sscal,sswap,snrm2
83: c fortran abs,amax1,min0,sqrt
84: c
85: c internal variables
86: c
87: integer j,jp,l,lp1,lup,maxj,pl,pu
88: real maxnrm,snrm2,tt
89: real sdot,nrmxl,t
90: logical negj,swapj
91: c
92: c
93: pl = 1
94: pu = 0
95: if (job .eq. 0) go to 60
96: c
97: c pivoting has been requested. rearrange the columns
98: c according to jpvt.
99: c
100: do 20 j = 1, p
101: swapj = jpvt(j) .gt. 0
102: negj = jpvt(j) .lt. 0
103: jpvt(j) = j
104: if (negj) jpvt(j) = -j
105: if (.not.swapj) go to 10
106: if (j .ne. pl) call sswap(n,x(1,pl),1,x(1,j),1)
107: jpvt(j) = jpvt(pl)
108: jpvt(pl) = j
109: pl = pl + 1
110: 10 continue
111: 20 continue
112: pu = p
113: do 50 jj = 1, p
114: j = p - jj + 1
115: if (jpvt(j) .ge. 0) go to 40
116: jpvt(j) = -jpvt(j)
117: if (j .eq. pu) go to 30
118: call sswap(n,x(1,pu),1,x(1,j),1)
119: jp = jpvt(pu)
120: jpvt(pu) = jpvt(j)
121: jpvt(j) = jp
122: 30 continue
123: pu = pu - 1
124: 40 continue
125: 50 continue
126: 60 continue
127: c
128: c compute the norms of the free columns.
129: c
130: if (pu .lt. pl) go to 80
131: do 70 j = pl, pu
132: qraux(j) = snrm2(n,x(1,j),1)
133: work(j) = qraux(j)
134: 70 continue
135: 80 continue
136: c
137: c perform the householder reduction of x.
138: c
139: lup = min0(n,p)
140: do 200 l = 1, lup
141: if (l .lt. pl .or. l .ge. pu) go to 120
142: c
143: c locate the column of largest norm and bring it
144: c into the pivot position.
145: c
146: maxnrm = 0.0e0
147: maxj = l
148: do 100 j = l, pu
149: if (qraux(j) .le. maxnrm) go to 90
150: maxnrm = qraux(j)
151: maxj = j
152: 90 continue
153: 100 continue
154: if (maxj .eq. l) go to 110
155: call sswap(n,x(1,l),1,x(1,maxj),1)
156: qraux(maxj) = qraux(l)
157: work(maxj) = work(l)
158: jp = jpvt(maxj)
159: jpvt(maxj) = jpvt(l)
160: jpvt(l) = jp
161: 110 continue
162: 120 continue
163: qraux(l) = 0.0e0
164: if (l .eq. n) go to 190
165: c
166: c compute the householder transformation for column l.
167: c
168: nrmxl = snrm2(n-l+1,x(l,l),1)
169: if (nrmxl .eq. 0.0e0) go to 180
170: if (x(l,l) .ne. 0.0e0) nrmxl = sign(nrmxl,x(l,l))
171: call sscal(n-l+1,1.0e0/nrmxl,x(l,l),1)
172: x(l,l) = 1.0e0 + x(l,l)
173: c
174: c apply the transformation to the remaining columns,
175: c updating the norms.
176: c
177: lp1 = l + 1
178: if (p .lt. lp1) go to 170
179: do 160 j = lp1, p
180: t = -sdot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
181: call saxpy(n-l+1,t,x(l,l),1,x(l,j),1)
182: if (j .lt. pl .or. j .gt. pu) go to 150
183: if (qraux(j) .eq. 0.0e0) go to 150
184: tt = 1.0e0 - (abs(x(l,j))/qraux(j))**2
185: tt = amax1(tt,0.0e0)
186: t = tt
187: tt = 1.0e0 + 0.05e0*tt*(qraux(j)/work(j))**2
188: if (tt .eq. 1.0e0) go to 130
189: qraux(j) = qraux(j)*sqrt(t)
190: go to 140
191: 130 continue
192: qraux(j) = snrm2(n-l,x(l+1,j),1)
193: work(j) = qraux(j)
194: 140 continue
195: 150 continue
196: 160 continue
197: 170 continue
198: c
199: c save the transformation.
200: c
201: qraux(l) = x(l,l)
202: x(l,l) = -nrmxl
203: 180 continue
204: 190 continue
205: 200 continue
206: return
207: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.