89
90
91
92
93
94
95 INTEGER LDA, LWORK, M, N
96
97
98 DOUBLE PRECISION A( LDA, * ), S( * ), WORK( LWORK )
99
100
101
102
103
104 DOUBLE PRECISION ZERO, ONE
105 parameter( zero = 0.0d0, one = 1.0d0 )
106
107
108 INTEGER I, INFO, ISCL, J, MN
109 DOUBLE PRECISION ANRM, BIGNUM, NRMSVL, SMLNUM
110
111
112 DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DNRM2
114
115
117
118
119 INTRINSIC dble, max, min
120
121
122 DOUBLE PRECISION DUMMY( 1 )
123
124
125
127
128
129
130 IF( lwork.LT.max( m*n+4*min( m, n )+max( m, n ),
131 $ m*n+2*min( m, n )+4*n) ) THEN
132 CALL xerbla(
'DQRT12', 7 )
133 RETURN
134 END IF
135
136
137
138 mn = min( m, n )
139 IF( mn.LE.zero )
140 $ RETURN
141
142 nrmsvl =
dnrm2( mn, s, 1 )
143
144
145
146 CALL dlaset(
'Full', m, n, zero, zero, work, m )
147 DO j = 1, n
148 DO i = 1, min( j, m )
149 work( ( j-1 )*m+i ) = a( i, j )
150 END DO
151 END DO
152
153
154
156 bignum = one / smlnum
157
158
159
160 anrm =
dlange(
'M', m, n, work, m, dummy )
161 iscl = 0
162 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
163
164
165
166 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, work, m, info )
167 iscl = 1
168 ELSE IF( anrm.GT.bignum ) THEN
169
170
171
172 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, work, m, info )
173 iscl = 1
174 END IF
175
176 IF( anrm.NE.zero ) THEN
177
178
179
180 CALL dgebd2( m, n, work, m, work( m*n+1 ), work( m*n+mn+1 ),
181 $ work( m*n+2*mn+1 ), work( m*n+3*mn+1 ),
182 $ work( m*n+4*mn+1 ), info )
183 CALL dbdsqr(
'Upper', mn, 0, 0, 0, work( m*n+1 ),
184 $ work( m*n+mn+1 ), dummy, mn, dummy, 1, dummy, mn,
185 $ work( m*n+2*mn+1 ), info )
186
187 IF( iscl.EQ.1 ) THEN
188 IF( anrm.GT.bignum ) THEN
189 CALL dlascl(
'G', 0, 0, bignum, anrm, mn, 1,
190 $ work( m*n+1 ), mn, info )
191 END IF
192 IF( anrm.LT.smlnum ) THEN
193 CALL dlascl(
'G', 0, 0, smlnum, anrm, mn, 1,
194 $ work( m*n+1 ), mn, info )
195 END IF
196 END IF
197
198 ELSE
199
200 DO i = 1, mn
201 work( m*n+i ) = zero
202 END DO
203 END IF
204
205
206
207 CALL daxpy( mn, -one, s, 1, work( m*n+1 ), 1 )
208
210 $ (
dlamch(
'Epsilon') * dble( max( m, n ) ) )
211
212 IF( nrmsvl.NE.zero )
214
215 RETURN
216
217
218
subroutine xerbla(srname, info)
double precision function dqrt12(m, n, a, lda, s, work, lwork)
DQRT12
double precision function dasum(n, dx, incx)
DASUM
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
subroutine dbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DBDSQR
subroutine dgebd2(m, n, a, lda, d, e, tauq, taup, work, info)
DGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
double precision function dlamch(cmach)
DLAMCH
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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 dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
real(wp) function dnrm2(n, x, incx)
DNRM2