152
153
154
155
156
157
158 LOGICAL UPPER
159 REAL A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
160 $ SNU, SNV
161
162
163
164
165
166 REAL ZERO
167 parameter( zero = 0.0e+0 )
168
169
170 REAL A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
171 $ AVB21, AVB22, CSL, CSR, D, S1, S2, SNL,
172 $ SNR, UA11R, UA22R, VB11R, VB22R, B, C, R, UA11,
173 $ UA12, UA21, UA22, VB11, VB12, VB21, VB22
174
175
177
178
179 INTRINSIC abs
180
181
182
183 IF( upper ) THEN
184
185
186
187
188
189
190 a = a1*b3
191 d = a3*b1
192 b = a2*b1 - a1*b2
193
194
195
196
197
198
199 CALL slasv2( a, b, d, s1, s2, snr, csr, snl, csl )
200
201 IF( abs( csl ).GE.abs( snl ) .OR. abs( csr ).GE.abs( snr ) )
202 $ THEN
203
204
205
206
207 ua11r = csl*a1
208 ua12 = csl*a2 + snl*a3
209
210 vb11r = csr*b1
211 vb12 = csr*b2 + snr*b3
212
213 aua12 = abs( csl )*abs( a2 ) + abs( snl )*abs( a3 )
214 avb12 = abs( csr )*abs( b2 ) + abs( snr )*abs( b3 )
215
216
217
218 IF( ( abs( ua11r )+abs( ua12 ) ).NE.zero ) THEN
219 IF( aua12 / ( abs( ua11r )+abs( ua12 ) ).LE.avb12 /
220 $ ( abs( vb11r )+abs( vb12 ) ) ) THEN
221 CALL slartg( -ua11r, ua12, csq, snq, r )
222 ELSE
223 CALL slartg( -vb11r, vb12, csq, snq, r )
224 END IF
225 ELSE
226 CALL slartg( -vb11r, vb12, csq, snq, r )
227 END IF
228
229 csu = csl
230 snu = -snl
231 csv = csr
232 snv = -snr
233
234 ELSE
235
236
237
238
239 ua21 = -snl*a1
240 ua22 = -snl*a2 + csl*a3
241
242 vb21 = -snr*b1
243 vb22 = -snr*b2 + csr*b3
244
245 aua22 = abs( snl )*abs( a2 ) + abs( csl )*abs( a3 )
246 avb22 = abs( snr )*abs( b2 ) + abs( csr )*abs( b3 )
247
248
249
250 IF( ( abs( ua21 )+abs( ua22 ) ).NE.zero ) THEN
251 IF( aua22 / ( abs( ua21 )+abs( ua22 ) ).LE.avb22 /
252 $ ( abs( vb21 )+abs( vb22 ) ) ) THEN
253 CALL slartg( -ua21, ua22, csq, snq, r )
254 ELSE
255 CALL slartg( -vb21, vb22, csq, snq, r )
256 END IF
257 ELSE
258 CALL slartg( -vb21, vb22, csq, snq, r )
259 END IF
260
261 csu = snl
262 snu = csl
263 csv = snr
264 snv = csr
265
266 END IF
267
268 ELSE
269
270
271
272
273
274
275 a = a1*b3
276 d = a3*b1
277 c = a2*b3 - a3*b2
278
279
280
281
282
283
284 CALL slasv2( a, c, d, s1, s2, snr, csr, snl, csl )
285
286 IF( abs( csr ).GE.abs( snr ) .OR. abs( csl ).GE.abs( snl ) )
287 $ THEN
288
289
290
291
292 ua21 = -snr*a1 + csr*a2
293 ua22r = csr*a3
294
295 vb21 = -snl*b1 + csl*b2
296 vb22r = csl*b3
297
298 aua21 = abs( snr )*abs( a1 ) + abs( csr )*abs( a2 )
299 avb21 = abs( snl )*abs( b1 ) + abs( csl )*abs( b2 )
300
301
302
303 IF( ( abs( ua21 )+abs( ua22r ) ).NE.zero ) THEN
304 IF( aua21 / ( abs( ua21 )+abs( ua22r ) ).LE.avb21 /
305 $ ( abs( vb21 )+abs( vb22r ) ) ) THEN
306 CALL slartg( ua22r, ua21, csq, snq, r )
307 ELSE
308 CALL slartg( vb22r, vb21, csq, snq, r )
309 END IF
310 ELSE
311 CALL slartg( vb22r, vb21, csq, snq, r )
312 END IF
313
314 csu = csr
315 snu = -snr
316 csv = csl
317 snv = -snl
318
319 ELSE
320
321
322
323
324 ua11 = csr*a1 + snr*a2
325 ua12 = snr*a3
326
327 vb11 = csl*b1 + snl*b2
328 vb12 = snl*b3
329
330 aua11 = abs( csr )*abs( a1 ) + abs( snr )*abs( a2 )
331 avb11 = abs( csl )*abs( b1 ) + abs( snl )*abs( b2 )
332
333
334
335 IF( ( abs( ua11 )+abs( ua12 ) ).NE.zero ) THEN
336 IF( aua11 / ( abs( ua11 )+abs( ua12 ) ).LE.avb11 /
337 $ ( abs( vb11 )+abs( vb12 ) ) ) THEN
338 CALL slartg( ua12, ua11, csq, snq, r )
339 ELSE
340 CALL slartg( vb12, vb11, csq, snq, r )
341 END IF
342 ELSE
343 CALL slartg( vb12, vb11, csq, snq, r )
344 END IF
345
346 csu = snr
347 snu = csr
348 csv = snl
349 snv = csl
350
351 END IF
352
353 END IF
354
355 RETURN
356
357
358
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine slasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.