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