LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dgelsd()

subroutine dgelsd ( integer  M,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( * )  S,
double precision  RCOND,
integer  RANK,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DGELSD computes the minimum-norm solution to a linear least squares problem for GE matrices

Download DGELSD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DGELSD computes the minimum-norm solution to a real linear least
 squares problem:
     minimize 2-norm(| b - A*x |)
 using the singular value decomposition (SVD) of A. A is an M-by-N
 matrix which may be rank-deficient.

 Several right hand side vectors b and solution vectors x can be
 handled in a single call; they are stored as the columns of the
 M-by-NRHS right hand side matrix B and the N-by-NRHS solution
 matrix X.

 The problem is solved in three steps:
 (1) Reduce the coefficient matrix A to bidiagonal form with
     Householder transformations, reducing the original problem
     into a "bidiagonal least squares problem" (BLS)
 (2) Solve the BLS using a divide and conquer approach.
 (3) Apply back all the Householder transformations to solve
     the original least squares problem.

 The effective rank of A is determined by treating as zero those
 singular values which are less than RCOND times the largest singular
 value.

 The divide and conquer algorithm makes very mild assumptions about
 floating point arithmetic. It will work on machines with a guard
 digit in add/subtract, or on those binary machines without guard
 digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
 Cray-2. It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.
Parameters
[in]M
          M is INTEGER
          The number of rows of A. M >= 0.
[in]N
          N is INTEGER
          The number of columns of A. N >= 0.
[in]NRHS
          NRHS is INTEGER
          The number of right hand sides, i.e., the number of columns
          of the matrices B and X. NRHS >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit, A has been destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
          On entry, the M-by-NRHS right hand side matrix B.
          On exit, B is overwritten by the N-by-NRHS solution
          matrix X.  If m >= n and RANK = n, the residual
          sum-of-squares for the solution in the i-th column is given
          by the sum of squares of elements n+1:m in that column.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,max(M,N)).
[out]S
          S is DOUBLE PRECISION array, dimension (min(M,N))
          The singular values of A in decreasing order.
          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
[in]RCOND
          RCOND is DOUBLE PRECISION
          RCOND is used to determine the effective rank of A.
          Singular values S(i) <= RCOND*S(1) are treated as zero.
          If RCOND < 0, machine precision is used instead.
[out]RANK
          RANK is INTEGER
          The effective rank of A, i.e., the number of singular values
          which are greater than RCOND*S(1).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK must be at least 1.
          The exact minimum amount of workspace needed depends on M,
          N and NRHS. As long as LWORK is at least
              12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
          if M is greater than or equal to N or
              12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
          if M is less than N, the code will execute correctly.
          SMLSIZ is returned by ILAENV and is equal to the maximum
          size of the subproblems at the bottom of the computation
          tree (usually about 25), and
             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          LIWORK >= max(1, 3 * MINMN * NLVL + 11 * MINMN),
          where MINMN = MIN( M,N ).
          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  the algorithm for computing the SVD failed to converge;
                if INFO = i, i off-diagonal elements of an intermediate
                bidiagonal form did not converge to zero.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 211 of file dgelsd.f.

211 *
212 * -- LAPACK driver routine (version 3.7.1) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215 * June 2017
216 *
217 * .. Scalar Arguments ..
218  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
219  DOUBLE PRECISION rcond
220 * ..
221 * .. Array Arguments ..
222  INTEGER iwork( * )
223  DOUBLE PRECISION a( lda, * ), b( ldb, * ), s( * ), work( * )
224 * ..
225 *
226 * =====================================================================
227 *
228 * .. Parameters ..
229  DOUBLE PRECISION zero, one, two
230  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
231 * ..
232 * .. Local Scalars ..
233  LOGICAL lquery
234  INTEGER iascl, ibscl, ie, il, itau, itaup, itauq,
235  $ ldwork, liwork, maxmn, maxwrk, minmn, minwrk,
236  $ mm, mnthr, nlvl, nwork, smlsiz, wlalsd
237  DOUBLE PRECISION anrm, bignum, bnrm, eps, sfmin, smlnum
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dgebrd, dgelqf, dgeqrf, dlabad, dlacpy, dlalsd,
242 * ..
243 * .. External Functions ..
244  INTEGER ilaenv
245  DOUBLE PRECISION dlamch, dlange
246  EXTERNAL ilaenv, dlamch, dlange
247 * ..
248 * .. Intrinsic Functions ..
249  INTRINSIC dble, int, log, max, min
250 * ..
251 * .. Executable Statements ..
252 *
253 * Test the input arguments.
254 *
255  info = 0
256  minmn = min( m, n )
257  maxmn = max( m, n )
258  mnthr = ilaenv( 6, 'DGELSD', ' ', m, n, nrhs, -1 )
259  lquery = ( lwork.EQ.-1 )
260  IF( m.LT.0 ) THEN
261  info = -1
262  ELSE IF( n.LT.0 ) THEN
263  info = -2
264  ELSE IF( nrhs.LT.0 ) THEN
265  info = -3
266  ELSE IF( lda.LT.max( 1, m ) ) THEN
267  info = -5
268  ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
269  info = -7
270  END IF
271 *
272  smlsiz = ilaenv( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
273 *
274 * Compute workspace.
275 * (Note: Comments in the code beginning "Workspace:" describe the
276 * minimal amount of workspace needed at that point in the code,
277 * as well as the preferred amount for good performance.
278 * NB refers to the optimal block size for the immediately
279 * following subroutine, as returned by ILAENV.)
280 *
281  minwrk = 1
282  liwork = 1
283  minmn = max( 1, minmn )
284  nlvl = max( int( log( dble( minmn ) / dble( smlsiz+1 ) ) /
285  $ log( two ) ) + 1, 0 )
286 *
287  IF( info.EQ.0 ) THEN
288  maxwrk = 0
289  liwork = 3*minmn*nlvl + 11*minmn
290  mm = m
291  IF( m.GE.n .AND. m.GE.mnthr ) THEN
292 *
293 * Path 1a - overdetermined, with many more rows than columns.
294 *
295  mm = n
296  maxwrk = max( maxwrk, n+n*ilaenv( 1, 'DGEQRF', ' ', m, n,
297  $ -1, -1 ) )
298  maxwrk = max( maxwrk, n+nrhs*
299  $ ilaenv( 1, 'DORMQR', 'LT', m, nrhs, n, -1 ) )
300  END IF
301  IF( m.GE.n ) THEN
302 *
303 * Path 1 - overdetermined or exactly determined.
304 *
305  maxwrk = max( maxwrk, 3*n+( mm+n )*
306  $ ilaenv( 1, 'DGEBRD', ' ', mm, n, -1, -1 ) )
307  maxwrk = max( maxwrk, 3*n+nrhs*
308  $ ilaenv( 1, 'DORMBR', 'QLT', mm, nrhs, n, -1 ) )
309  maxwrk = max( maxwrk, 3*n+( n-1 )*
310  $ ilaenv( 1, 'DORMBR', 'PLN', n, nrhs, n, -1 ) )
311  wlalsd = 9*n+2*n*smlsiz+8*n*nlvl+n*nrhs+(smlsiz+1)**2
312  maxwrk = max( maxwrk, 3*n+wlalsd )
313  minwrk = max( 3*n+mm, 3*n+nrhs, 3*n+wlalsd )
314  END IF
315  IF( n.GT.m ) THEN
316  wlalsd = 9*m+2*m*smlsiz+8*m*nlvl+m*nrhs+(smlsiz+1)**2
317  IF( n.GE.mnthr ) THEN
318 *
319 * Path 2a - underdetermined, with many more columns
320 * than rows.
321 *
322  maxwrk = m + m*ilaenv( 1, 'DGELQF', ' ', m, n, -1, -1 )
323  maxwrk = max( maxwrk, m*m+4*m+2*m*
324  $ ilaenv( 1, 'DGEBRD', ' ', m, m, -1, -1 ) )
325  maxwrk = max( maxwrk, m*m+4*m+nrhs*
326  $ ilaenv( 1, 'DORMBR', 'QLT', m, nrhs, m, -1 ) )
327  maxwrk = max( maxwrk, m*m+4*m+( m-1 )*
328  $ ilaenv( 1, 'DORMBR', 'PLN', m, nrhs, m, -1 ) )
329  IF( nrhs.GT.1 ) THEN
330  maxwrk = max( maxwrk, m*m+m+m*nrhs )
331  ELSE
332  maxwrk = max( maxwrk, m*m+2*m )
333  END IF
334  maxwrk = max( maxwrk, m+nrhs*
335  $ ilaenv( 1, 'DORMLQ', 'LT', n, nrhs, m, -1 ) )
336  maxwrk = max( maxwrk, m*m+4*m+wlalsd )
337 ! XXX: Ensure the Path 2a case below is triggered. The workspace
338 ! calculation should use queries for all routines eventually.
339  maxwrk = max( maxwrk,
340  $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
341  ELSE
342 *
343 * Path 2 - remaining underdetermined cases.
344 *
345  maxwrk = 3*m + ( n+m )*ilaenv( 1, 'DGEBRD', ' ', m, n,
346  $ -1, -1 )
347  maxwrk = max( maxwrk, 3*m+nrhs*
348  $ ilaenv( 1, 'DORMBR', 'QLT', m, nrhs, n, -1 ) )
349  maxwrk = max( maxwrk, 3*m+m*
350  $ ilaenv( 1, 'DORMBR', 'PLN', n, nrhs, m, -1 ) )
351  maxwrk = max( maxwrk, 3*m+wlalsd )
352  END IF
353  minwrk = max( 3*m+nrhs, 3*m+m, 3*m+wlalsd )
354  END IF
355  minwrk = min( minwrk, maxwrk )
356  work( 1 ) = maxwrk
357  iwork( 1 ) = liwork
358 
359  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
360  info = -12
361  END IF
362  END IF
363 *
364  IF( info.NE.0 ) THEN
365  CALL xerbla( 'DGELSD', -info )
366  RETURN
367  ELSE IF( lquery ) THEN
368  GO TO 10
369  END IF
370 *
371 * Quick return if possible.
372 *
373  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
374  rank = 0
375  RETURN
376  END IF
377 *
378 * Get machine parameters.
379 *
380  eps = dlamch( 'P' )
381  sfmin = dlamch( 'S' )
382  smlnum = sfmin / eps
383  bignum = one / smlnum
384  CALL dlabad( smlnum, bignum )
385 *
386 * Scale A if max entry outside range [SMLNUM,BIGNUM].
387 *
388  anrm = dlange( 'M', m, n, a, lda, work )
389  iascl = 0
390  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
391 *
392 * Scale matrix norm up to SMLNUM.
393 *
394  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
395  iascl = 1
396  ELSE IF( anrm.GT.bignum ) THEN
397 *
398 * Scale matrix norm down to BIGNUM.
399 *
400  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
401  iascl = 2
402  ELSE IF( anrm.EQ.zero ) THEN
403 *
404 * Matrix all zero. Return zero solution.
405 *
406  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
407  CALL dlaset( 'F', minmn, 1, zero, zero, s, 1 )
408  rank = 0
409  GO TO 10
410  END IF
411 *
412 * Scale B if max entry outside range [SMLNUM,BIGNUM].
413 *
414  bnrm = dlange( 'M', m, nrhs, b, ldb, work )
415  ibscl = 0
416  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
417 *
418 * Scale matrix norm up to SMLNUM.
419 *
420  CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
421  ibscl = 1
422  ELSE IF( bnrm.GT.bignum ) THEN
423 *
424 * Scale matrix norm down to BIGNUM.
425 *
426  CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
427  ibscl = 2
428  END IF
429 *
430 * If M < N make sure certain entries of B are zero.
431 *
432  IF( m.LT.n )
433  $ CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
434 *
435 * Overdetermined case.
436 *
437  IF( m.GE.n ) THEN
438 *
439 * Path 1 - overdetermined or exactly determined.
440 *
441  mm = m
442  IF( m.GE.mnthr ) THEN
443 *
444 * Path 1a - overdetermined, with many more rows than columns.
445 *
446  mm = n
447  itau = 1
448  nwork = itau + n
449 *
450 * Compute A=Q*R.
451 * (Workspace: need 2*N, prefer N+N*NB)
452 *
453  CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
454  $ lwork-nwork+1, info )
455 *
456 * Multiply B by transpose(Q).
457 * (Workspace: need N+NRHS, prefer N+NRHS*NB)
458 *
459  CALL dormqr( 'L', 'T', m, nrhs, n, a, lda, work( itau ), b,
460  $ ldb, work( nwork ), lwork-nwork+1, info )
461 *
462 * Zero out below R.
463 *
464  IF( n.GT.1 ) THEN
465  CALL dlaset( 'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
466  END IF
467  END IF
468 *
469  ie = 1
470  itauq = ie + n
471  itaup = itauq + n
472  nwork = itaup + n
473 *
474 * Bidiagonalize R in A.
475 * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
476 *
477  CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
478  $ work( itaup ), work( nwork ), lwork-nwork+1,
479  $ info )
480 *
481 * Multiply B by transpose of left bidiagonalizing vectors of R.
482 * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
483 *
484  CALL dormbr( 'Q', 'L', 'T', mm, nrhs, n, a, lda, work( itauq ),
485  $ b, ldb, work( nwork ), lwork-nwork+1, info )
486 *
487 * Solve the bidiagonal least squares problem.
488 *
489  CALL dlalsd( 'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
490  $ rcond, rank, work( nwork ), iwork, info )
491  IF( info.NE.0 ) THEN
492  GO TO 10
493  END IF
494 *
495 * Multiply B by right bidiagonalizing vectors of R.
496 *
497  CALL dormbr( 'P', 'L', 'N', n, nrhs, n, a, lda, work( itaup ),
498  $ b, ldb, work( nwork ), lwork-nwork+1, info )
499 *
500  ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
501  $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
502 *
503 * Path 2a - underdetermined, with many more columns than rows
504 * and sufficient workspace for an efficient algorithm.
505 *
506  ldwork = m
507  IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
508  $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
509  itau = 1
510  nwork = m + 1
511 *
512 * Compute A=L*Q.
513 * (Workspace: need 2*M, prefer M+M*NB)
514 *
515  CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
516  $ lwork-nwork+1, info )
517  il = nwork
518 *
519 * Copy L to WORK(IL), zeroing out above its diagonal.
520 *
521  CALL dlacpy( 'L', m, m, a, lda, work( il ), ldwork )
522  CALL dlaset( 'U', m-1, m-1, zero, zero, work( il+ldwork ),
523  $ ldwork )
524  ie = il + ldwork*m
525  itauq = ie + m
526  itaup = itauq + m
527  nwork = itaup + m
528 *
529 * Bidiagonalize L in WORK(IL).
530 * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
531 *
532  CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
533  $ work( itauq ), work( itaup ), work( nwork ),
534  $ lwork-nwork+1, info )
535 *
536 * Multiply B by transpose of left bidiagonalizing vectors of L.
537 * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
538 *
539  CALL dormbr( 'Q', 'L', 'T', m, nrhs, m, work( il ), ldwork,
540  $ work( itauq ), b, ldb, work( nwork ),
541  $ lwork-nwork+1, info )
542 *
543 * Solve the bidiagonal least squares problem.
544 *
545  CALL dlalsd( 'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
546  $ rcond, rank, work( nwork ), iwork, info )
547  IF( info.NE.0 ) THEN
548  GO TO 10
549  END IF
550 *
551 * Multiply B by right bidiagonalizing vectors of L.
552 *
553  CALL dormbr( 'P', 'L', 'N', m, nrhs, m, work( il ), ldwork,
554  $ work( itaup ), b, ldb, work( nwork ),
555  $ lwork-nwork+1, info )
556 *
557 * Zero out below first M rows of B.
558 *
559  CALL dlaset( 'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
560  nwork = itau + m
561 *
562 * Multiply transpose(Q) by B.
563 * (Workspace: need M+NRHS, prefer M+NRHS*NB)
564 *
565  CALL dormlq( 'L', 'T', n, nrhs, m, a, lda, work( itau ), b,
566  $ ldb, work( nwork ), lwork-nwork+1, info )
567 *
568  ELSE
569 *
570 * Path 2 - remaining underdetermined cases.
571 *
572  ie = 1
573  itauq = ie + m
574  itaup = itauq + m
575  nwork = itaup + m
576 *
577 * Bidiagonalize A.
578 * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
579 *
580  CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
581  $ work( itaup ), work( nwork ), lwork-nwork+1,
582  $ info )
583 *
584 * Multiply B by transpose of left bidiagonalizing vectors.
585 * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
586 *
587  CALL dormbr( 'Q', 'L', 'T', m, nrhs, n, a, lda, work( itauq ),
588  $ b, ldb, work( nwork ), lwork-nwork+1, info )
589 *
590 * Solve the bidiagonal least squares problem.
591 *
592  CALL dlalsd( 'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
593  $ rcond, rank, work( nwork ), iwork, info )
594  IF( info.NE.0 ) THEN
595  GO TO 10
596  END IF
597 *
598 * Multiply B by right bidiagonalizing vectors of A.
599 *
600  CALL dormbr( 'P', 'L', 'N', n, nrhs, m, a, lda, work( itaup ),
601  $ b, ldb, work( nwork ), lwork-nwork+1, info )
602 *
603  END IF
604 *
605 * Undo scaling.
606 *
607  IF( iascl.EQ.1 ) THEN
608  CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
609  CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
610  $ info )
611  ELSE IF( iascl.EQ.2 ) THEN
612  CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
613  CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
614  $ info )
615  END IF
616  IF( ibscl.EQ.1 ) THEN
617  CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
618  ELSE IF( ibscl.EQ.2 ) THEN
619  CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
620  END IF
621 *
622  10 CONTINUE
623  work( 1 ) = maxwrk
624  iwork( 1 ) = liwork
625  RETURN
626 *
627 * End of DGELSD
628 *
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
DLALSD uses the singular value decomposition of A to solve the least squares problem.
Definition: dlalsd.f:181
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
Definition: dgebrd.f:207
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.
Definition: dlascl.f:145
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...
Definition: dlaset.f:112
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
Definition: dgelqf.f:137
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:169
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 ...
Definition: dlange.f:116
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR
Definition: dormbr.f:197
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
Definition: dormlq.f:169
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:138
Here is the call graph for this function:
Here is the caller graph for this function: