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
118
119
120 INTRINSIC dble, max, min
121
122
123 DOUBLE PRECISION DUMMY( 1 )
124
125
126
128
129
130
131 IF( lwork.LT.max( m*n+4*min( m, n )+max( m, n ),
132 $ m*n+2*min( m, n )+4*n) ) THEN
133 CALL xerbla(
'DQRT12', 7 )
134 RETURN
135 END IF
136
137
138
139 mn = min( m, n )
140 IF( mn.LE.zero )
141 $ RETURN
142
143 nrmsvl =
dnrm2( mn, s, 1 )
144
145
146
147 CALL dlaset(
'Full', m, n, zero, zero, work, m )
148 DO 20 j = 1, n
149 DO 10 i = 1, min( j, m )
150 work( ( j-1 )*m+i ) = a( i, j )
151 10 CONTINUE
152 20 CONTINUE
153
154
155
157 bignum = one / smlnum
158 CALL dlabad( smlnum, bignum )
159
160
161
162 anrm =
dlange(
'M', m, n, work, m, dummy )
163 iscl = 0
164 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
165
166
167
168 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, work, m, info )
169 iscl = 1
170 ELSE IF( anrm.GT.bignum ) THEN
171
172
173
174 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, work, m, info )
175 iscl = 1
176 END IF
177
178 IF( anrm.NE.zero ) THEN
179
180
181
182 CALL dgebd2( m, n, work, m, work( m*n+1 ), work( m*n+mn+1 ),
183 $ work( m*n+2*mn+1 ), work( m*n+3*mn+1 ),
184 $ work( m*n+4*mn+1 ), info )
185 CALL dbdsqr(
'Upper', mn, 0, 0, 0, work( m*n+1 ),
186 $ work( m*n+mn+1 ), dummy, mn, dummy, 1, dummy, mn,
187 $ work( m*n+2*mn+1 ), info )
188
189 IF( iscl.EQ.1 ) THEN
190 IF( anrm.GT.bignum ) THEN
191 CALL dlascl(
'G', 0, 0, bignum, anrm, mn, 1,
192 $ work( m*n+1 ), mn, info )
193 END IF
194 IF( anrm.LT.smlnum ) THEN
195 CALL dlascl(
'G', 0, 0, smlnum, anrm, mn, 1,
196 $ work( m*n+1 ), mn, info )
197 END IF
198 END IF
199
200 ELSE
201
202 DO 30 i = 1, mn
203 work( m*n+i ) = zero
204 30 CONTINUE
205 END IF
206
207
208
209 CALL daxpy( mn, -one, s, 1, work( m*n+1 ), 1 )
211 $ (
dlamch(
'Epsilon' )*dble( max( m, n ) ) )
212 IF( nrmsvl.NE.zero )
214
215 RETURN
216
217
218
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 dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
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 dasum(N, DX, INCX)
DASUM
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
double precision function dqrt12(M, N, A, LDA, S, WORK, LWORK)
DQRT12
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 dgebd2(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
DGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
real(wp) function dnrm2(n, x, incx)
DNRM2