116
117
118
119
120
121
122 CHARACTER TRANS
123 INTEGER LDA, LDX, LWORK, M, N, NRHS
124
125
126 REAL A( LDA, * ), WORK( LWORK ), X( LDX, * )
127
128
129
130
131
132 REAL ZERO, ONE
133 parameter( zero = 0.0e0, one = 1.0e0 )
134
135
136 LOGICAL TPSD
137 INTEGER I, INFO, J, LDWORK
138 REAL ANRM, ERR, XNRM
139
140
141 REAL RWORK( 1 )
142
143
144 LOGICAL LSAME
145 REAL SLAMCH, SLANGE
147
148
150
151
152 INTRINSIC abs, max, min, real
153
154
155
157 IF(
lsame( trans,
'N' ) )
THEN
158 ldwork = m + nrhs
159 tpsd = .false.
160 IF( lwork.LT.( m+nrhs )*( n+2 ) ) THEN
161 CALL xerbla(
'SQRT14', 10 )
162 RETURN
163 ELSE IF( n.LE.0 .OR. nrhs.LE.0 ) THEN
164 RETURN
165 END IF
166 ELSE IF(
lsame( trans,
'T' ) )
THEN
167 ldwork = m
168 tpsd = .true.
169 IF( lwork.LT.( n+nrhs )*( m+2 ) ) THEN
170 CALL xerbla(
'SQRT14', 10 )
171 RETURN
172 ELSE IF( m.LE.0 .OR. nrhs.LE.0 ) THEN
173 RETURN
174 END IF
175 ELSE
176 CALL xerbla(
'SQRT14', 1 )
177 RETURN
178 END IF
179
180
181
182 CALL slacpy(
'All', m, n, a, lda, work, ldwork )
183 anrm =
slange(
'M', m, n, work, ldwork, rwork )
184 IF( anrm.NE.zero )
185 $
CALL slascl(
'G', 0, 0, anrm, one, m, n, work, ldwork, info )
186
187
188
189 IF( tpsd ) THEN
190
191
192
193 CALL slacpy(
'All', m, nrhs, x, ldx, work( n*ldwork+1 ),
194 $ ldwork )
195 xnrm =
slange(
'M', m, nrhs, work( n*ldwork+1 ), ldwork,
196 $ rwork )
197 IF( xnrm.NE.zero )
198 $
CALL slascl(
'G', 0, 0, xnrm, one, m, nrhs,
199 $ work( n*ldwork+1 ), ldwork, info )
200
201
202
203 CALL sgeqr2( m, n+nrhs, work, ldwork,
204 $ work( ldwork*( n+nrhs )+1 ),
205 $ work( ldwork*( n+nrhs )+min( m, n+nrhs )+1 ),
206 $ info )
207
208
209
210
211 err = zero
212 DO 20 j = n + 1, n + nrhs
213 DO 10 i = n + 1, min( m, j )
214 err = max( err, abs( work( i+( j-1 )*m ) ) )
215 10 CONTINUE
216 20 CONTINUE
217
218 ELSE
219
220
221
222 DO 40 i = 1, n
223 DO 30 j = 1, nrhs
224 work( m+j+( i-1 )*ldwork ) = x( i, j )
225 30 CONTINUE
226 40 CONTINUE
227
228 xnrm =
slange(
'M', nrhs, n, work( m+1 ), ldwork, rwork )
229 IF( xnrm.NE.zero )
230 $
CALL slascl(
'G', 0, 0, xnrm, one, nrhs, n, work( m+1 ),
231 $ ldwork, info )
232
233
234
235 CALL sgelq2( ldwork, n, work, ldwork, work( ldwork*n+1 ),
236 $ work( ldwork*( n+1 )+1 ), info )
237
238
239
240
241 err = zero
242 DO 60 j = m + 1, n
243 DO 50 i = j, ldwork
244 err = max( err, abs( work( i+( j-1 )*ldwork ) ) )
245 50 CONTINUE
246 60 CONTINUE
247
248 END IF
249
250 sqrt14 = err / ( real( max( m, n, nrhs ) )*
slamch(
'Epsilon' ) )
251
252 RETURN
253
254
255
subroutine xerbla(srname, info)
subroutine sgelq2(m, n, a, lda, tau, work, info)
SGELQ2 computes the LQ factorization of a general rectangular matrix using an unblocked algorithm.
subroutine sgeqr2(m, n, a, lda, tau, work, info)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
logical function lsame(ca, cb)
LSAME
real function sqrt14(trans, m, n, nrhs, a, lda, x, ldx, work, lwork)
SQRT14