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