97
98
99
100
101
102
103 INTEGER LDA, LWORK, M, N
104
105
106 DOUBLE PRECISION RWORK( * ), S( * )
107 COMPLEX*16 A( LDA, * ), WORK( LWORK )
108
109
110
111
112
113 DOUBLE PRECISION ZERO, ONE
114 parameter( zero = 0.0d0, one = 1.0d0 )
115
116
117 INTEGER I, INFO, ISCL, J, MN
118 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM
119
120
121 DOUBLE PRECISION DUMMY( 1 )
122
123
124 DOUBLE PRECISION DASUM, DLAMCH, DNRM2, ZLANGE
126
127
130
131
132 INTRINSIC dble, dcmplx, max, min
133
134
135
137
138
139
140 IF( lwork.LT.m*n+2*min( m, n )+max( m, n ) ) THEN
141 CALL xerbla(
'ZQRT12', 7 )
142 RETURN
143 END IF
144
145
146
147 mn = min( m, n )
148 IF( mn.LE.zero )
149 $ RETURN
150
151 nrmsvl =
dnrm2( mn, s, 1 )
152
153
154
155 CALL zlaset(
'Full', m, n, dcmplx( zero ), dcmplx( zero ), work,
156 $ m )
157 DO 20 j = 1, n
158 DO 10 i = 1, min( j, m )
159 work( ( j-1 )*m+i ) = a( i, j )
160 10 CONTINUE
161 20 CONTINUE
162
163
164
166 bignum = one / smlnum
167 CALL dlabad( smlnum, bignum )
168
169
170
171 anrm =
zlange(
'M', m, n, work, m, dummy )
172 iscl = 0
173 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
174
175
176
177 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, work, m, info )
178 iscl = 1
179 ELSE IF( anrm.GT.bignum ) THEN
180
181
182
183 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, work, m, info )
184 iscl = 1
185 END IF
186
187 IF( anrm.NE.zero ) THEN
188
189
190
191 CALL zgebd2( m, n, work, m, rwork( 1 ), rwork( mn+1 ),
192 $ work( m*n+1 ), work( m*n+mn+1 ),
193 $ work( m*n+2*mn+1 ), info )
194 CALL dbdsqr(
'Upper', mn, 0, 0, 0, rwork( 1 ), rwork( mn+1 ),
195 $ dummy, mn, dummy, 1, dummy, mn, rwork( 2*mn+1 ),
196 $ info )
197
198 IF( iscl.EQ.1 ) THEN
199 IF( anrm.GT.bignum ) THEN
200 CALL dlascl(
'G', 0, 0, bignum, anrm, mn, 1, rwork( 1 ),
201 $ mn, info )
202 END IF
203 IF( anrm.LT.smlnum ) THEN
204 CALL dlascl(
'G', 0, 0, smlnum, anrm, mn, 1, rwork( 1 ),
205 $ mn, info )
206 END IF
207 END IF
208
209 ELSE
210
211 DO 30 i = 1, mn
212 rwork( i ) = zero
213 30 CONTINUE
214 END IF
215
216
217
218 CALL daxpy( mn, -one, s, 1, rwork( 1 ), 1 )
220 $ (
dlamch(
'Epsilon' )*dble( max( m, n ) ) )
221 IF( nrmsvl.NE.zero )
223
224 RETURN
225
226
227
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
double precision function zqrt12(M, N, A, LDA, S, WORK, LWORK, RWORK)
ZQRT12
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zgebd2(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
ZGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
double precision function dasum(N, DX, INCX)
DASUM
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
real(wp) function dnrm2(n, x, incx)
DNRM2