1 /*
2 * $Id: dbdsqr.c,v 1.1 2005-09-18 22:04:39 dhmunro Exp $
3 * LAPACK matrix solver using SVD.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 #include "dg.h"
12
13
14 /*---blas routines---*/
15 /* dscal dswap drot */
16
17
18
19 /*---converted nutty string switches to single characters (lower case)---*/
20 #define lsame(x,y) ((x)==(y))
21
22
23
24 /*-----Fortran intrinsics converted-----*/
25 extern double pow(double,double); /* used only to take 1/8th power */
26 #define abs(x) ((x)>=0?(x):-(x))
27 extern double sqrt(double);
28 #define dble(x) ((double)x)
29 #define sign(x,y) ((((x)<0)!=((y)<0))?-(x):(x))
30 #define min(x,y) ((x)<(y)?(x):(y))
31 #define max(x,y) ((x)>(y)?(x):(y))
32 /*-----end of Fortran intrinsics-----*/
33
34
dbdsqr(char uplo,long n,long ncvt,long nru,long ncc,double d[],double e[],double vt[],long ldvt,double u[],long ldu,double c[],long ldc,double work[],long * info)35 void dbdsqr( char uplo, long n, long ncvt, long nru, long ncc,
36 double d[], double e[], double vt[], long ldvt,
37 double u[],long ldu, double c[], long ldc,
38 double work[], long *info )
39 {
40 /**
41 * -- LAPACK routine (version 1.1) --
42 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
43 * Courant Institute, Argonne National Lab, and Rice University
44 * March 31, 1993
45 *
46 * .. Scalar Arguments ..*/
47 /** ..
48 * .. Array Arguments ..*/
49 #undef work_1
50 #define work_1(a1) work[a1-1]
51 #undef vt_2
52 #define vt_2(a1,a2) vt[a1-1+ldvt*(a2-1)]
53 #undef u_2
54 #define u_2(a1,a2) u[a1-1+ldu*(a2-1)]
55 #undef e_1
56 #define e_1(a1) e[a1-1]
57 #undef d_1
58 #define d_1(a1) d[a1-1]
59 #undef c_2
60 #define c_2(a1,a2) c[a1-1+ldc*(a2-1)]
61 /** ..
62 *
63 * Purpose
64 * =======
65 *
66 * DBDSQR computes the singular value decomposition (SVD) of a real
67 * N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P'
68 * denotes the transpose of P), where S is a diagonal matrix with
69 * non-negative diagonal elements (the singular values of B), and Q
70 * and P are orthogonal matrices.
71 *
72 * The routine computes S, and optionally computes U * Q, P' * VT,
73 * or Q' * C, for given real input matrices U, VT, and C.
74 *
75 * See "Computing Small Singular Values of Bidiagonal Matrices With
76 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
77 * LAPACK Working Note #3, for a detailed description of the algorithm.
78 *
79 * Arguments
80 * =========
81 *
82 * UPLO (input) CHARACTER*1
83 * = 'U': B is upper bidiagonal;
84 * = 'L': B is lower bidiagonal.
85 *
86 * N (input) INTEGER
87 * The order of the matrix B. N >= 0.
88 *
89 * NCVT (input) INTEGER
90 * The number of columns of the matrix VT. NCVT >= 0.
91 *
92 * NRU (input) INTEGER
93 * The number of rows of the matrix U. NRU >= 0.
94 *
95 * NCC (input) INTEGER
96 * The number of columns of the matrix C. NCC >= 0.
97 *
98 * D (input/output) DOUBLE PRECISION array, dimension (N)
99 * On entry, the n diagonal elements of the bidiagonal matrix B.
100 * On exit, if INFO=0, the singular values of B in decreasing
101 * order.
102 *
103 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
104 * On entry, the (n-1) off-diagonal elements of the bidiagonal
105 * matrix B.
106 * On normal exit, E is destroyed.
107 *
108 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
109 * On entry, an N-by-NCVT matrix VT.
110 * On exit, VT is overwritten by P' * VT.
111 * VT is not referenced if NCVT = 0.
112 *
113 * LDVT (input) INTEGER
114 * The leading dimension of the array VT.
115 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
116 *
117 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
118 * On entry, an NRU-by-N matrix U.
119 * On exit, U is overwritten by U * Q.
120 * U is not referenced if NRU = 0.
121 *
122 * LDU (input) INTEGER
123 * The leading dimension of the array U. LDU >= max(1,NRU).
124 *
125 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
126 * On entry, an N-by-NCC matrix C.
127 * On exit, C is overwritten by Q' * C.
128 * C is not referenced if NCC = 0.
129 *
130 * LDC (input) INTEGER
131 * The leading dimension of the array C.
132 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
133 *
134 * WORK (workspace) DOUBLE PRECISION array, dimension
135 * (MAX( 1, 4*N-4 ))
136 * WORK is not referenced if NCVT = NRU = NCC = 0.
137 *
138 * INFO (output) INTEGER
139 * = 0: successful exit
140 * < 0: If INFO = -i, the i-th argument had an illegal value
141 * > 0: the algorithm did not converge; D and E contain the
142 * elements of a bidiagonal matrix which is orthogonally
143 * similar to the input matrix B; if INFO = i, i
144 * elements of E have not converged to zero.
145 *
146 * Internal Parameters
147 * ===================
148 *
149 * TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
150 * TOLMUL controls the convergence criterion of the QR loop.
151 * If it is positive, TOLMUL*EPS is the desired relative
152 * precision in the computed singular values.
153 * If it is negative, abs(TOLMUL*EPS*sigma_max) is the
154 * desired absolute accuracy in the computed singular
155 * values (corresponds to relative accuracy
156 * abs(TOLMUL*EPS) in the largest singular value.
157 * abs(TOLMUL) should be between 1 and 1/EPS, and preferably
158 * between 10 (for fast convergence) and .1/EPS
159 * (for there to be some accuracy in the results).
160 * Default is to lose at either one eighth or 2 of the
161 * available decimal digits in each computed singular value
162 * (whichever is smaller).
163 *
164 * MAXITR INTEGER, default = 6
165 * MAXITR controls the maximum number of passes of the
166 * algorithm through its inner loop. The algorithms stops
167 * (and so fails to converge) if the number of passes
168 * through the inner loop exceeds MAXITR*N**2.
169 *
170 * =====================================================================
171 *
172 * .. Parameters ..*/
173 #undef zero
174 #define zero 0.0e0
175 #undef one
176 #define one 1.0e0
177 #undef negone
178 #define negone -1.0e0
179 #undef hndrth
180 #define hndrth 0.01e0
181 #undef ten
182 #define ten 10.0e0
183 #undef hndrd
184 #define hndrd 100.0e0
185 #undef meigth
186 #define meigth -0.125e0
187 #undef maxitr
188 #define maxitr 6
189 /** ..
190 * .. Local Scalars ..*/
191 int rotate;
192 long i, idir=0, irot, isub, iter, iuplo, j, job, ll=0,
193 lll, m, maxit, nm1, nm12, nm13, oldll, oldm;
194 double abse, abss, cosl, cosr, cs, eps, f, g, gap,
195 gmax, h, mu, oldcs, oldsn, r, shift, sigmn,
196 sigmx, sinl, sinr, sll, smax, smin, sminl,
197 sminlo=0.0, sminoa=0.0, sn, thresh, tol, tolmul, unfl;
198 /** ..
199 * .. Intrinsic Functions ..*/
200 /* intrinsic abs, dble, max, min, sign, sqrt;*/
201 /** ..
202 * .. Executable Statements ..
203 *
204 * Test the input parameters.
205 **/
206 /*-----implicit-declarations-----*/
207 /*-----end-of-declarations-----*/
208 *info = 0;
209 iuplo = 0;
210 if( lsame( uplo, 'u' ) )
211 iuplo = 1;
212 if( lsame( uplo, 'l' ) )
213 iuplo = 2;
214 if( iuplo==0 ) {
215 *info = -1;
216 } else if( n<0 ) {
217 *info = -2;
218 } else if( ncvt<0 ) {
219 *info = -3;
220 } else if( nru<0 ) {
221 *info = -4;
222 } else if( ncc<0 ) {
223 *info = -5;
224 } else if( ( ncvt==0 && ldvt<1 ) ||
225 ( ncvt>0 && ldvt<max( 1, n ) ) ) {
226 *info = -9;
227 } else if( ldu<max( 1, nru ) ) {
228 *info = -11;
229 } else if( ( ncc==0 && ldc<1 ) ||
230 ( ncc>0 && ldc<max( 1, n ) ) ) {
231 *info = -13;
232 }
233 if( *info!=0 ) {
234 xerbla( "dbdsqr", -*info );
235 return;
236 }
237 if( n==0 )
238 return;
239 if( n==1 )
240 goto L_190;
241 /**
242 * ROTATE is true if any singular vectors desired, false otherwise
243 **/
244 rotate = ( ncvt>0 ) || ( nru>0 ) || ( ncc>0 );
245 nm1 = n - 1;
246 nm12 = nm1 + nm1;
247 nm13 = nm12 + nm1;
248 /**
249 * Get machine constants
250 **/
251 eps = dlamch( 'e'/*psilon*/ );
252 unfl = dlamch( 's'/*afe minimum*/ );
253 tolmul = max( ten, min( hndrd, pow(eps,meigth) ) );
254 tol = tolmul*eps;
255 /**
256 * If matrix lower bidiagonal, rotate to be upper bidiagonal
257 * by applying Givens rotations on the left
258 **/
259 if( iuplo==2 ) {
260 for (i=1 ; i<=n - 1 ; i+=1) {
261 dlartg( d_1( i ), e_1( i ), &cs, &sn, &r );
262 d_1( i ) = r;
263 e_1( i ) = sn*d_1( i+1 );
264 d_1( i+1 ) = cs*d_1( i+1 );
265 if( rotate ) {
266 work_1( i ) = cs;
267 work_1( nm1+i ) = sn;
268 }
269 }
270 /**
271 * Update singular vectors if desired
272 **/
273 if( nru>0 )
274 dlasr( 'r', 'v', 'f', nru, n, &work_1( 1 ), &work_1( n ), u,
275 ldu );
276 if( ncc>0 )
277 dlasr( 'l', 'v', 'f', n, ncc, &work_1( 1 ), &work_1( n ), c,
278 ldc );
279 }
280 /**
281 * Compute approximate maximum, minimum singular values
282 **/
283 smax = abs( d_1( n ) );
284 for (i=1 ; i<=n - 1 ; i+=1) {
285 smax = max( smax, max(abs( d_1( i ) ), abs( e_1( i ) )) );
286 }
287 sminl = zero;
288 if( tol>=zero ) {
289 sminoa = abs( d_1( 1 ) );
290 if( sminoa==zero )
291 goto L_40;
292 mu = sminoa;
293 for (i=2 ; i<=n ; i+=1) {
294 mu = abs( d_1( i ) )*( mu / ( mu+abs( e_1( i-1 ) ) ) );
295 sminoa = min( sminoa, mu );
296 if( sminoa==zero )
297 goto L_40;
298 }
299 L_40:
300 sminoa = sminoa / sqrt( dble( n ) );
301 }
302 /**
303 * Prepare for main iteration loop for the singular values
304 **/
305 maxit = maxitr*n*n;
306 iter = 0;
307 oldll = -1;
308 oldm = -1;
309 if( ncc==0 && nru==0 && ncvt==0 ) {
310 /**
311 * No singular vectors desired
312 **/
313 job = 0;
314 } else {
315 /**
316 * Singular vectors desired
317 **/
318 job = 1;
319 }
320 if( tol>=zero ) {
321 /**
322 * Relative accuracy desired
323 **/
324 thresh = max( tol*sminoa, maxit*unfl );
325 } else {
326 /**
327 * Absolute accuracy desired
328 **/
329 thresh = max( abs( tol )*smax, maxit*unfl );
330 }
331 /**
332 * M points to last entry of unconverged part of matrix
333 **/
334 m = n;
335 /**
336 * Begin main iteration loop
337 **/
338 L_50:
339 /**
340 * Check for convergence or exceeding iteration count
341 **/
342 if( m<=1 )
343 goto L_190;
344 if( iter>maxit )
345 goto L_230;
346 /**
347 * Find diagonal block of matrix to work on
348 **/
349 if( tol<zero && abs( d_1( m ) )<=thresh )
350 d_1( m ) = zero;
351 smax = abs( d_1( m ) );
352 smin = smax;
353 for (lll=1 ; lll<=m ; lll+=1) {
354 ll = m - lll;
355 if( ll==0 )
356 goto L_80;
357 abss = abs( d_1( ll ) );
358 abse = abs( e_1( ll ) );
359 if( tol<zero && abss<=thresh )
360 d_1( ll ) = zero;
361 if( abse<=thresh )
362 goto L_70;
363 smin = min( smin, abss );
364 smax = max( smax, max(abss, abse) );
365 }
366 L_70:
367 e_1( ll ) = zero;
368 /**
369 * Matrix splits since E(LL) = 0
370 **/
371 if( ll==m-1 ) {
372 /**
373 * Convergence of bottom singular value, return to top of loop
374 **/
375 m = m - 1;
376 goto L_50;
377 }
378 L_80:
379 ll = ll + 1;
380 /**
381 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
382 **/
383 if( ll==m-1 ) {
384 /**
385 * 2 by 2 block, handle separately
386 **/
387 dlasv2( d_1( m-1 ), e_1( m-1 ), d_1( m ), &sigmn, &sigmx, &sinr,
388 &cosr, &sinl, &cosl );
389 d_1( m-1 ) = sigmx;
390 e_1( m-1 ) = zero;
391 d_1( m ) = sigmn;
392 /**
393 * Compute singular vectors, if desired
394 **/
395 if( ncvt>0 )
396 cblas_drot( ncvt, &vt_2( m-1, 1 ), ldvt, &vt_2( m, 1 ), ldvt, cosr,
397 sinr );
398 if( nru>0 )
399 cblas_drot( nru, &u_2( 1, m-1 ), 1, &u_2( 1, m ), 1, cosl, sinl );
400 if( ncc>0 )
401 cblas_drot( ncc, &c_2( m-1, 1 ), ldc, &c_2( m, 1 ), ldc, cosl,
402 sinl );
403 m = m - 2;
404 goto L_50;
405 }
406 /**
407 * If working on new submatrix, choose shift direction
408 * (from larger end diagonal entry towards smaller)
409 **/
410 if( ll>oldm || m<oldll ) {
411 if( abs( d_1( ll ) )>=abs( d_1( m ) ) ) {
412 /**
413 * Chase bulge from top (big end) to bottom (small end)
414 **/
415 idir = 1;
416 } else {
417 /**
418 * Chase bulge from bottom (big end) to top (small end)
419 **/
420 idir = 2;
421 }
422 }
423 /**
424 * Apply convergence tests
425 **/
426 if( idir==1 ) {
427 /**
428 * Run convergence test in forward direction
429 * First apply standard test to bottom of matrix
430 **/
431 if( abs( e_1( m-1 ) )<=abs( tol )*abs( d_1( m ) ) ||
432 ( tol<zero && abs( e_1( m-1 ) )<=thresh ) ) {
433 e_1( m-1 ) = zero;
434 goto L_50;
435 }
436
437 if( tol>=zero ) {
438 /**
439 * If relative accuracy desired,
440 * apply convergence criterion forward
441 **/
442 mu = abs( d_1( ll ) );
443 sminl = mu;
444 for (lll=ll ; lll<=m - 1 ; lll+=1) {
445 if( abs( e_1( lll ) )<=tol*mu ) {
446 e_1( lll ) = zero;
447 goto L_50;
448 }
449 sminlo = sminl;
450 mu = abs( d_1( lll+1 ) )*( mu / ( mu+abs( e_1( lll ) ) ) );
451 sminl = min( sminl, mu );
452 }
453 /**
454 * If singular values only wanted, apply gap test to bottom
455 * end of matrix
456 **/
457 if( job==0 ) {
458 gap = sminlo / sqrt( dble( m-ll ) ) - abs( d_1( m ) );
459 if( gap>zero ) {
460 abss = abs( d_1( m ) );
461 abse = abs( e_1( m-1 ) );
462 gmax = max( gap, max(abss, abse) );
463 if( ( abse / gmax )*( abse / gmax )<=tol*( gap / gmax )*
464 ( abss / gmax ) ) {
465 e_1( m-1 ) = zero;
466 goto L_50;
467 }
468 }
469 }
470 }
471 } else {
472 /**
473 * Run convergence test in backward direction
474 * First apply standard test to top of matrix
475 **/
476 if( abs( e_1( ll ) )<=abs( tol )*abs( d_1( ll ) ) ||
477 ( tol<zero && abs( e_1( ll ) )<=thresh ) ) {
478 e_1( ll ) = zero;
479 goto L_50;
480 }
481
482 if( tol>=zero ) {
483 /**
484 * If relative accuracy desired,
485 * apply convergence criterion backward
486 **/
487 mu = abs( d_1( m ) );
488 sminl = mu;
489 for (lll=m - 1 ; lll>=ll ; lll+=-1) {
490 if( abs( e_1( lll ) )<=tol*mu ) {
491 e_1( lll ) = zero;
492 goto L_50;
493 }
494 sminlo = sminl;
495 mu = abs( d_1( lll ) )*( mu / ( mu+abs( e_1( lll ) ) ) );
496 sminl = min( sminl, mu );
497 }
498 /**
499 * If singular values only wanted, apply gap test to top
500 * end of matrix
501 **/
502 if( job==0 ) {
503 gap = sminlo / sqrt( dble( m-ll ) ) - abs( d_1( ll ) );
504 if( gap>zero ) {
505 abss = abs( d_1( ll ) );
506 abse = abs( e_1( ll ) );
507 gmax = max( gap, max(abss, abse) );
508 if( ( abse / gmax )*( abse / gmax )<=tol*( gap / gmax )*
509 ( abss / gmax ) ) {
510 e_1( ll ) = zero;
511 goto L_50;
512 }
513 }
514 }
515 }
516 }
517 oldll = ll;
518 oldm = m;
519 /**
520 * Compute shift. First, test if shifting would ruin relative
521 * accuracy, and if so set the shift to zero.
522 **/
523 if( tol>=zero && n*tol*( sminl / smax )<=
524 max( eps, hndrth*tol ) ) {
525 /**
526 * Use a zero shift to avoid loss of relative accuracy
527 **/
528 shift = zero;
529 } else {
530 /**
531 * Compute the shift from 2-by-2 block at end of matrix
532 **/
533 if( idir==1 ) {
534 sll = abs( d_1( ll ) );
535 dlas2( d_1( m-1 ), e_1( m-1 ), d_1( m ), &shift, &r );
536 } else {
537 sll = abs( d_1( m ) );
538 dlas2( d_1( ll ), e_1( ll ), d_1( ll+1 ), &shift, &r );
539 }
540 /**
541 * Test if shift negligible, and if so set to zero
542 **/
543 if( sll>zero ) {
544 if( ( shift / sll )*( shift / sll )<eps )
545 shift = zero;
546 }
547 }
548 /**
549 * Increment iteration count
550 **/
551 iter = iter + m - ll;
552 /**
553 * If SHIFT = 0, do simplified QR iteration
554 **/
555 if( shift==zero ) {
556 if( idir==1 ) {
557 /**
558 * Chase bulge from top to bottom
559 **/
560 cs = one;
561 oldcs = one;
562 /**
563 * Save cosines and sines if singular vectors desired
564 **/
565 if( rotate ) {
566
567 dlartg( d_1( ll )*cs, e_1( ll ), &cs, &sn, &r );
568 dlartg( oldcs*r, d_1( ll+1 )*sn, &oldcs, &oldsn,
569 &d_1( ll ) );
570 work_1( 1 ) = cs;
571 work_1( 1+nm1 ) = sn;
572 work_1( 1+nm12 ) = oldcs;
573 work_1( 1+nm13 ) = oldsn;
574 irot = 1;
575 for (i=ll + 1 ; i<=m - 1 ; i+=1) {
576 dlartg( d_1( i )*cs, e_1( i ), &cs, &sn, &r );
577 e_1( i-1 ) = oldsn*r;
578 dlartg( oldcs*r, d_1( i+1 )*sn, &oldcs, &oldsn,
579 &d_1( i ) );
580 irot = irot + 1;
581 work_1( irot ) = cs;
582 work_1( irot+nm1 ) = sn;
583 work_1( irot+nm12 ) = oldcs;
584 work_1( irot+nm13 ) = oldsn;
585 }
586 h = d_1( m )*cs;
587 d_1( m ) = h*oldcs;
588 e_1( m-1 ) = h*oldsn;
589 /**
590 * Update singular vectors
591 **/
592 if( ncvt>0 )
593 dlasr( 'l', 'v', 'f', m-ll+1, ncvt, &work_1( 1 ),
594 &work_1( n ), &vt_2( ll, 1 ), ldvt );
595 if( nru>0 )
596 dlasr( 'r', 'v', 'f', nru, m-ll+1,
597 &work_1( nm12+1 ), &work_1( nm13+1 ),
598 &u_2( 1, ll ), ldu );
599 if( ncc>0 )
600 dlasr( 'l', 'v', 'f', m-ll+1, ncc,
601 &work_1( nm12+1 ), &work_1( nm13+1 ),
602 &c_2( ll, 1 ), ldc );
603
604 } else {
605
606 dlartg( d_1( ll )*cs, e_1( ll ), &cs, &sn, &r );
607 dlartg( oldcs*r, d_1( ll+1 )*sn, &oldcs, &oldsn,
608 &d_1( ll ) );
609 for (i=ll + 1 ; i<=m - 1 ; i+=1) {
610 dlartg( d_1( i )*cs, e_1( i ), &cs, &sn, &r );
611 e_1( i-1 ) = oldsn*r;
612 dlartg( oldcs*r, d_1( i+1 )*sn, &oldcs, &oldsn,
613 &d_1( i ) );
614 }
615 h = d_1( m )*cs;
616 d_1( m ) = h*oldcs;
617 e_1( m-1 ) = h*oldsn;
618
619 }
620 /**
621 * Test convergence
622 **/
623 if( abs( e_1( m-1 ) )<=thresh )
624 e_1( m-1 ) = zero;
625
626 } else {
627 /**
628 * Chase bulge from bottom to top
629 **/
630 cs = one;
631 oldcs = one;
632 /**
633 * Save cosines and sines if singular vectors desired
634 **/
635 if( rotate ) {
636
637 dlartg( d_1( m )*cs, e_1( m-1 ), &cs, &sn, &r );
638 dlartg( oldcs*r, d_1( m-1 )*sn, &oldcs, &oldsn, &d_1( m ) );
639 work_1( m-ll ) = cs;
640 work_1( m-ll+nm1 ) = -sn;
641 work_1( m-ll+nm12 ) = oldcs;
642 work_1( m-ll+nm13 ) = -oldsn;
643 irot = m - ll;
644 for (i=m - 1 ; i>=ll + 1 ; i+=-1) {
645 dlartg( d_1( i )*cs, e_1( i-1 ), &cs, &sn, &r );
646 e_1( i ) = oldsn*r;
647 dlartg( oldcs*r, d_1( i-1 )*sn, &oldcs, &oldsn,
648 &d_1( i ) );
649 irot = irot - 1;
650 work_1( irot ) = cs;
651 work_1( irot+nm1 ) = -sn;
652 work_1( irot+nm12 ) = oldcs;
653 work_1( irot+nm13 ) = -oldsn;
654 }
655 h = d_1( ll )*cs;
656 d_1( ll ) = h*oldcs;
657 e_1( ll ) = h*oldsn;
658 /**
659 * Update singular vectors
660 **/
661 if( ncvt>0 )
662 dlasr( 'l', 'v', 'b', m-ll+1, ncvt,
663 &work_1( nm12+1 ), &work_1( nm13+1 ),
664 &vt_2( ll, 1 ), ldvt );
665 if( nru>0 )
666 dlasr( 'r', 'v', 'b', nru, m-ll+1, &work_1( 1 ),
667 &work_1( n ), &u_2( 1, ll ), ldu );
668 if( ncc>0 )
669 dlasr( 'l', 'v', 'b', m-ll+1, ncc, &work_1( 1 ),
670 &work_1( n ), &c_2( ll, 1 ), ldc );
671
672 } else {
673
674 dlartg( d_1( m )*cs, e_1( m-1 ), &cs, &sn, &r );
675 dlartg( oldcs*r, d_1( m-1 )*sn, &oldcs, &oldsn, &d_1( m ) );
676 for (i=m - 1 ; i>=ll + 1 ; i+=-1) {
677 dlartg( d_1( i )*cs, e_1( i-1 ), &cs, &sn, &r );
678 e_1( i ) = oldsn*r;
679 dlartg( oldcs*r, d_1( i-1 )*sn, &oldcs, &oldsn,
680 &d_1( i ) );
681 }
682 h = d_1( ll )*cs;
683 d_1( ll ) = h*oldcs;
684 e_1( ll ) = h*oldsn;
685
686 }
687 /**
688 * Test convergence
689 **/
690 if( abs( e_1( ll ) )<=thresh )
691 e_1( ll ) = zero;
692 }
693 } else {
694 /**
695 * Use nonzero shift
696 **/
697 if( idir==1 ) {
698 /**
699 * Chase bulge from top to bottom
700 **/
701 f = ( abs( d_1( ll ) )-shift )*
702 ( sign( one, d_1( ll ) )+shift / d_1( ll ) );
703 g = e_1( ll );
704 /**
705 * Save cosines and sines if singular vectors desired
706 **/
707 if( rotate ) {
708
709 dlartg( f, g, &cosr, &sinr, &r );
710 f = cosr*d_1( ll ) + sinr*e_1( ll );
711 e_1( ll ) = cosr*e_1( ll ) - sinr*d_1( ll );
712 g = sinr*d_1( ll+1 );
713 d_1( ll+1 ) = cosr*d_1( ll+1 );
714 dlartg( f, g, &cosl, &sinl, &r );
715 d_1( ll ) = r;
716 f = cosl*e_1( ll ) + sinl*d_1( ll+1 );
717 d_1( ll+1 ) = cosl*d_1( ll+1 ) - sinl*e_1( ll );
718 g = sinl*e_1( ll+1 );
719 e_1( ll+1 ) = cosl*e_1( ll+1 );
720 work_1( 1 ) = cosr;
721 work_1( 1+nm1 ) = sinr;
722 work_1( 1+nm12 ) = cosl;
723 work_1( 1+nm13 ) = sinl;
724 irot = 1;
725 for (i=ll + 1 ; i<=m - 2 ; i+=1) {
726 dlartg( f, g, &cosr, &sinr, &r );
727 e_1( i-1 ) = r;
728 f = cosr*d_1( i ) + sinr*e_1( i );
729 e_1( i ) = cosr*e_1( i ) - sinr*d_1( i );
730 g = sinr*d_1( i+1 );
731 d_1( i+1 ) = cosr*d_1( i+1 );
732 dlartg( f, g, &cosl, &sinl, &r );
733 d_1( i ) = r;
734 f = cosl*e_1( i ) + sinl*d_1( i+1 );
735 d_1( i+1 ) = cosl*d_1( i+1 ) - sinl*e_1( i );
736 g = sinl*e_1( i+1 );
737 e_1( i+1 ) = cosl*e_1( i+1 );
738 irot = irot + 1;
739 work_1( irot ) = cosr;
740 work_1( irot+nm1 ) = sinr;
741 work_1( irot+nm12 ) = cosl;
742 work_1( irot+nm13 ) = sinl;
743 }
744 dlartg( f, g, &cosr, &sinr, &r );
745 e_1( m-2 ) = r;
746 f = cosr*d_1( m-1 ) + sinr*e_1( m-1 );
747 e_1( m-1 ) = cosr*e_1( m-1 ) - sinr*d_1( m-1 );
748 g = sinr*d_1( m );
749 d_1( m ) = cosr*d_1( m );
750 dlartg( f, g, &cosl, &sinl, &r );
751 d_1( m-1 ) = r;
752 f = cosl*e_1( m-1 ) + sinl*d_1( m );
753 d_1( m ) = cosl*d_1( m ) - sinl*e_1( m-1 );
754 irot = irot + 1;
755 work_1( irot ) = cosr;
756 work_1( irot+nm1 ) = sinr;
757 work_1( irot+nm12 ) = cosl;
758 work_1( irot+nm13 ) = sinl;
759 e_1( m-1 ) = f;
760 /**
761 * Update singular vectors
762 **/
763 if( ncvt>0 )
764 dlasr( 'l', 'v', 'f', m-ll+1, ncvt, &work_1( 1 ),
765 &work_1( n ), &vt_2( ll, 1 ), ldvt );
766 if( nru>0 )
767 dlasr( 'r', 'v', 'f', nru, m-ll+1,
768 &work_1( nm12+1 ), &work_1( nm13+1 ),
769 &u_2( 1, ll ), ldu );
770 if( ncc>0 )
771 dlasr( 'l', 'v', 'f', m-ll+1, ncc,
772 &work_1( nm12+1 ), &work_1( nm13+1 ),
773 &c_2( ll, 1 ), ldc );
774
775 } else {
776
777 dlartg( f, g, &cosr, &sinr, &r );
778 f = cosr*d_1( ll ) + sinr*e_1( ll );
779 e_1( ll ) = cosr*e_1( ll ) - sinr*d_1( ll );
780 g = sinr*d_1( ll+1 );
781 d_1( ll+1 ) = cosr*d_1( ll+1 );
782 dlartg( f, g, &cosl, &sinl, &r );
783 d_1( ll ) = r;
784 f = cosl*e_1( ll ) + sinl*d_1( ll+1 );
785 d_1( ll+1 ) = cosl*d_1( ll+1 ) - sinl*e_1( ll );
786 g = sinl*e_1( ll+1 );
787 e_1( ll+1 ) = cosl*e_1( ll+1 );
788 for (i=ll + 1 ; i<=m - 2 ; i+=1) {
789 dlartg( f, g, &cosr, &sinr, &r );
790 e_1( i-1 ) = r;
791 f = cosr*d_1( i ) + sinr*e_1( i );
792 e_1( i ) = cosr*e_1( i ) - sinr*d_1( i );
793 g = sinr*d_1( i+1 );
794 d_1( i+1 ) = cosr*d_1( i+1 );
795 dlartg( f, g, &cosl, &sinl, &r );
796 d_1( i ) = r;
797 f = cosl*e_1( i ) + sinl*d_1( i+1 );
798 d_1( i+1 ) = cosl*d_1( i+1 ) - sinl*e_1( i );
799 g = sinl*e_1( i+1 );
800 e_1( i+1 ) = cosl*e_1( i+1 );
801 }
802 dlartg( f, g, &cosr, &sinr, &r );
803 e_1( m-2 ) = r;
804 f = cosr*d_1( m-1 ) + sinr*e_1( m-1 );
805 e_1( m-1 ) = cosr*e_1( m-1 ) - sinr*d_1( m-1 );
806 g = sinr*d_1( m );
807 d_1( m ) = cosr*d_1( m );
808 dlartg( f, g, &cosl, &sinl, &r );
809 d_1( m-1 ) = r;
810 f = cosl*e_1( m-1 ) + sinl*d_1( m );
811 d_1( m ) = cosl*d_1( m ) - sinl*e_1( m-1 );
812 e_1( m-1 ) = f;
813
814 }
815 /**
816 * Test convergence
817 **/
818 if( abs( e_1( m-1 ) )<=thresh )
819 e_1( m-1 ) = zero;
820
821 } else {
822 /**
823 * Chase bulge from bottom to top
824 **/
825 f = ( abs( d_1( m ) )-shift )*( sign( one, d_1( m ) )+shift /
826 d_1( m ) );
827 g = e_1( m-1 );
828 /**
829 * Save cosines and sines if singular vectors desired
830 **/
831 if( rotate ) {
832
833 dlartg( f, g, &cosr, &sinr, &r );
834 f = cosr*d_1( m ) + sinr*e_1( m-1 );
835 e_1( m-1 ) = cosr*e_1( m-1 ) - sinr*d_1( m );
836 g = sinr*d_1( m-1 );
837 d_1( m-1 ) = cosr*d_1( m-1 );
838 dlartg( f, g, &cosl, &sinl, &r );
839 d_1( m ) = r;
840 f = cosl*e_1( m-1 ) + sinl*d_1( m-1 );
841 d_1( m-1 ) = cosl*d_1( m-1 ) - sinl*e_1( m-1 );
842 g = sinl*e_1( m-2 );
843 e_1( m-2 ) = cosl*e_1( m-2 );
844 work_1( m-ll ) = cosr;
845 work_1( m-ll+nm1 ) = -sinr;
846 work_1( m-ll+nm12 ) = cosl;
847 work_1( m-ll+nm13 ) = -sinl;
848 irot = m - ll;
849 for (i=m - 1 ; i>=ll + 2 ; i+=-1) {
850 dlartg( f, g, &cosr, &sinr, &r );
851 e_1( i ) = r;
852 f = cosr*d_1( i ) + sinr*e_1( i-1 );
853 e_1( i-1 ) = cosr*e_1( i-1 ) - sinr*d_1( i );
854 g = sinr*d_1( i-1 );
855 d_1( i-1 ) = cosr*d_1( i-1 );
856 dlartg( f, g, &cosl, &sinl, &r );
857 d_1( i ) = r;
858 f = cosl*e_1( i-1 ) + sinl*d_1( i-1 );
859 d_1( i-1 ) = cosl*d_1( i-1 ) - sinl*e_1( i-1 );
860 g = sinl*e_1( i-2 );
861 e_1( i-2 ) = cosl*e_1( i-2 );
862 irot = irot - 1;
863 work_1( irot ) = cosr;
864 work_1( irot+nm1 ) = -sinr;
865 work_1( irot+nm12 ) = cosl;
866 work_1( irot+nm13 ) = -sinl;
867 }
868 dlartg( f, g, &cosr, &sinr, &r );
869 e_1( ll+1 ) = r;
870 f = cosr*d_1( ll+1 ) + sinr*e_1( ll );
871 e_1( ll ) = cosr*e_1( ll ) - sinr*d_1( ll+1 );
872 g = sinr*d_1( ll );
873 d_1( ll ) = cosr*d_1( ll );
874 dlartg( f, g, &cosl, &sinl, &r );
875 d_1( ll+1 ) = r;
876 f = cosl*e_1( ll ) + sinl*d_1( ll );
877 d_1( ll ) = cosl*d_1( ll ) - sinl*e_1( ll );
878 irot = irot - 1;
879 work_1( irot ) = cosr;
880 work_1( irot+nm1 ) = -sinr;
881 work_1( irot+nm12 ) = cosl;
882 work_1( irot+nm13 ) = -sinl;
883 e_1( ll ) = f;
884
885 } else {
886
887 dlartg( f, g, &cosr, &sinr, &r );
888 f = cosr*d_1( m ) + sinr*e_1( m-1 );
889 e_1( m-1 ) = cosr*e_1( m-1 ) - sinr*d_1( m );
890 g = sinr*d_1( m-1 );
891 d_1( m-1 ) = cosr*d_1( m-1 );
892 dlartg( f, g, &cosl, &sinl, &r );
893 d_1( m ) = r;
894 f = cosl*e_1( m-1 ) + sinl*d_1( m-1 );
895 d_1( m-1 ) = cosl*d_1( m-1 ) - sinl*e_1( m-1 );
896 g = sinl*e_1( m-2 );
897 e_1( m-2 ) = cosl*e_1( m-2 );
898 for (i=m - 1 ; i>=ll + 2 ; i+=-1) {
899 dlartg( f, g, &cosr, &sinr, &r );
900 e_1( i ) = r;
901 f = cosr*d_1( i ) + sinr*e_1( i-1 );
902 e_1( i-1 ) = cosr*e_1( i-1 ) - sinr*d_1( i );
903 g = sinr*d_1( i-1 );
904 d_1( i-1 ) = cosr*d_1( i-1 );
905 dlartg( f, g, &cosl, &sinl, &r );
906 d_1( i ) = r;
907 f = cosl*e_1( i-1 ) + sinl*d_1( i-1 );
908 d_1( i-1 ) = cosl*d_1( i-1 ) - sinl*e_1( i-1 );
909 g = sinl*e_1( i-2 );
910 e_1( i-2 ) = cosl*e_1( i-2 );
911 }
912 dlartg( f, g, &cosr, &sinr, &r );
913 e_1( ll+1 ) = r;
914 f = cosr*d_1( ll+1 ) + sinr*e_1( ll );
915 e_1( ll ) = cosr*e_1( ll ) - sinr*d_1( ll+1 );
916 g = sinr*d_1( ll );
917 d_1( ll ) = cosr*d_1( ll );
918 dlartg( f, g, &cosl, &sinl, &r );
919 d_1( ll+1 ) = r;
920 f = cosl*e_1( ll ) + sinl*d_1( ll );
921 d_1( ll ) = cosl*d_1( ll ) - sinl*e_1( ll );
922 e_1( ll ) = f;
923
924 }
925 /**
926 * Test convergence
927 **/
928 if( abs( e_1( ll ) )<=thresh )
929 e_1( ll ) = zero;
930 /**
931 * Update singular vectors if desired
932 **/
933 if( ncvt>0 )
934 dlasr( 'l', 'v', 'b', m-ll+1, ncvt, &work_1( nm12+1 ),
935 &work_1( nm13+1 ), &vt_2( ll, 1 ), ldvt );
936 if( nru>0 )
937 dlasr( 'r', 'v', 'b', nru, m-ll+1, &work_1( 1 ),
938 &work_1( n ), &u_2( 1, ll ), ldu );
939 if( ncc>0 )
940 dlasr( 'l', 'v', 'b', m-ll+1, ncc, &work_1( 1 ),
941 &work_1( n ), &c_2( ll, 1 ), ldc );
942 }
943 }
944 /**
945 * QR iteration finished, go back and check convergence
946 **/
947 goto L_50;
948 /**
949 * All singular values converged, so make them positive
950 **/
951 L_190:
952 for (i=1 ; i<=n ; i+=1) {
953 if( d_1( i )<zero ) {
954 d_1( i ) = -d_1( i );
955 /**
956 * Change sign of singular vectors, if desired
957 **/
958 if( ncvt>0 )
959 cblas_dscal( ncvt, negone, &vt_2( i, 1 ), ldvt );
960 }
961 }
962 /**
963 * Sort the singular values into decreasing order (insertion sort on
964 * singular values, but only one transposition per singular vector)
965 **/
966 for (i=1 ; i<=n - 1 ; i+=1) {
967 /**
968 * Scan for smallest D(I)
969 **/
970 isub = 1;
971 smin = d_1( 1 );
972 for (j=2 ; j<=n + 1 - i ; j+=1) {
973 if( d_1( j )<=smin ) {
974 isub = j;
975 smin = d_1( j );
976 }
977 }
978 if( isub!=n+1-i ) {
979 /**
980 * Swap singular values and vectors
981 **/
982 d_1( isub ) = d_1( n+1-i );
983 d_1( n+1-i ) = smin;
984 if( ncvt>0 )
985 cblas_dswap( ncvt, &vt_2( isub, 1 ), ldvt, &vt_2( n+1-i, 1 ),
986 ldvt );
987 if( nru>0 )
988 cblas_dswap( nru, &u_2( 1, isub ), 1, &u_2( 1, n+1-i ), 1 );
989 if( ncc>0 )
990 cblas_dswap( ncc, &c_2( isub, 1 ), ldc, &c_2( n+1-i, 1 ), ldc );
991 }
992 }
993 goto L_250;
994 /**
995 * Maximum number of iterations exceeded, failure to converge
996 **/
997 L_230:
998 *info = 0;
999 for (i=1 ; i<=n - 1 ; i+=1) {
1000 if( e_1( i )!=zero )
1001 *info = *info + 1;
1002 }
1003 L_250:
1004 return;
1005 /**
1006 * End of DBDSQR
1007 **/
1008 }
1009