scilab/modules/linear_algebra/src/fortran/intdgesv3.f

   1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 c Copyright (C) INRIA
   3 c$
   4 c This file must be used under the terms of the CeCILL.
   5 c This source file is licensed as described in the file COPYING, which
   6 c you should have received as part of this distribution.  The terms
   7 c are also available at
   8 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
   9 c$
  10       subroutine intdgesv3(fname)
  11
  12 c     a\b
  13
  14       include 'stack.h'
  15       logical getrhsvar,createvar
  16       logical checklhs,checkrhs
  17       character fname*(*)
  18       double precision ANORM, EPS, RCOND
  19       double precision dlamch, dlange
  20       integer vfinite
  21       external dlamch, dlange, vfinite
  22       intrinsic sqrt
  23 c
  24       minrhs=2
  25       maxrhs=2
  26       minlhs=1
  27       maxlhs=1
  28 c
  29       if(.not.checkrhs(fname,minrhs,maxrhs)) return
  30       if(.not.checklhs(fname,minlhs,maxlhs)) return
  31
  32       if(.not.getrhsvar(1,'d', MA, NA, lA)) return
  33       if(.not.getrhsvar(2,'d', MB, NRHS, lB)) return
  34       if(MB.eq.0) then
  35          lhsvar(1) = 2
  36          return
  37       endif
  38       if(MA .ne. MB) then
  39          call error(265)
  40          return
  41       endif
  42       M = MA
  43       N = NA
  44       if(M.eq.0 .or. N.eq.0) then
  45          if(.not.createvar(3,'d', 0, 0, lX)) return
  46          lhsvar(1) = 3
  47          return
  48       endif
  49 c     Check if A and B matrices contains Inf or NaN's
  50       if (vfinite(M*N, stk(lA)).eq.0) then
  51          call error(229)
  52          return
  53       endif
  54       if (vfinite(MB*NRHS, stk(lB)).eq.0) then
  55          call error(229)
  56          return
  57       endif
  58
  59
  60       if(.not.createvar(3,'d', M, N, lAF)) return
  61       if(.not.createvar(4,'d', N, NRHS, lX)) return
  62       if(.not.createvar(5,'d', max(M,N), NRHS, lXB)) return
  63       if(.not.createvar(6,'i', 1, 1, lRANK)) return
  64       if(.not.createvar(7,'i', 1, N, lIPIV)) return
  65       if(.not.createvar(8,'i', 1, N, lJPVT)) return
  66       if(.not.createvar(9,'i', 1, N, lIWORK)) return
  67       LWORKMIN = max( 4*N, min(M,N)+3*N+1, 2*min(M,N)+NRHS )
  68       LWORK=maxvol(10,'d')
  69       if(LWORK.le.LWORKMIN) then
  70          err=(LWORK-LWORKMIN)
  71          call error(17)
  72          return
  73       endif
  74       if(.not.createvar(10,'d',1,LWORK,lDWORK)) return
  75
  76       EPS = dlamch('eps')
  77       ANORM = dlange( '1', M, N, stk(lA), M, stk(lDWORK) )
  78 c
  79       if(M.eq.N) then
  80 c
  81 c     M = N
  82 c
  83          call DLACPY( 'F', N, N, stk(lA), N, stk(lAF), N )
  84 c     SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB )
  85          call DGETRF( N, N, stk(lAF), N, istk(lIPIV), INFO )
  86 c     SUBROUTINE DGETRF( N, N, A, LDA, IPIV, INFO )
  87          RCOND = 0.0d0
  88          if(INFO.eq.0) then
  89             call DGECON( '1', N, stk(lAF), N, ANORM, RCOND, stk(lDWORK),
  90      $           istk(lIWORK), INFO )
  91 c     SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK,
  92 c     $                        IWORK, INFO )
  93             if(RCOND.gt.sqrt(EPS)) then
  94                call DGETRS( 'N', N, NRHS, stk(lAF), N, istk(lIPIV),
  95      $              stk(lB), N, INFO )
  96 c     SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV,
  97 c     B, LDB, INFO )
  98                call DLACPY( 'F', N, NRHS, stk(lB), N, stk(lX), N )
  99 c     SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB )
 100                lhsvar(1) = 4
 101                return
 102             endif
 103          endif
 104          call writebufdgesv3(buf,RCOND)
 105          call msgs(5,1)
 106       endif
 107 c
 108 c     M.ne.N or A singular
 109 c
 110       RCOND = sqrt(EPS)
 111       call DLACPY( 'F', M, NRHS, stk(lB), M, stk(lXB), max(M,N) )
 112 c     SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB )
 113       do 10 i = 1, N
 114          istk(lJPVT+i-1) = 0
 115  10   continue
 116       call DGELSY1( M, N, NRHS, stk(lA), M, stk(lXB), max(M,N),
 117      $     istk(lJPVT), RCOND, istk(lRANK), stk(lDWORK),
 118      $     LWORK, INFO )
 119 c     SUBROUTINE DGELSY1( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
 120 c     $                    RANK, WORK, LWORK, INFO )
 121       if(info.ne.0) then
 122          return
 123       endif
 124       if( M.ne.N .and. istk(lRANK).lt.min(M,N) )then
 125          call msgs(9,istk(lRANK))
 126       endif
 127       call DLACPY( 'F', N, NRHS, stk(lXB), max(M,N), stk(lX), N )
 128 c     SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB )
 129
 130       lhsvar(1)=4
 131       return
 132 c
 133       end
 134
 135