scilab/modules/linear_algebra/src/fortran/intzgesv3.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
  11                 subroutine intzgesv3(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, zlange
  21       integer vfinite
  22       external dlamch, zlange, vfinite
  23       intrinsic sqrt
  24 c
  25       minrhs=2
  26       maxrhs=2
  27       minlhs=1
  28       maxlhs=1
  29 c
  30       if(.not.checkrhs(fname,minrhs,maxrhs)) return
  31       if(.not.checklhs(fname,minlhs,maxlhs)) return
  32
  33       if(.not.getrhsvar(1,'z', MA, NA, lA)) return
  34       if(.not.getrhsvar(2,'z', MB, NRHS, lB)) return
  35       if(MB.eq.0) then
  36          lhsvar(1) = 2
  37          return
  38       endif
  39       if(MA .ne. MB) then
  40          call error(265)
  41          return
  42       endif
  43       M = MA
  44       N = NA
  45       if(M.eq.0 .or. N.eq.0) then
  46          if(.not.createvar(3,'z', 0, 0, lX)) return
  47          lhsvar(1) = 3
  48          return
  49       endif
  50 c     Check if A and B matrices contains Inf or NaN's
  51       if (vfinite(M*N*2, zstk(lA)).eq.0) then
  52          call error(229)
  53          return
  54       endif
  55       if (vfinite(MB*NRHS*2, zstk(lB)).eq.0) then
  56          call error(229)
  57          return
  58       endif
  59
  60       if(.not.createvar(3,'z', M, N, lAF)) return
  61       if(.not.createvar(4,'z', N, NRHS, lX)) return
  62       if(.not.createvar(5,'z', 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,'d',1,2*N,lRWORK)) return
  67       LWORKMIN = max( 2*N, min(M,N)+max( 2*min(M,N), N+1,
  68      $     min(M,N)+NRHS ) )
  69       LWORK=maxvol(10,'z')
  70       if(LWORK.le.LWORKMIN) then
  71          err=2*(LWORK-LWORKMIN)
  72          call error(17)
  73          return
  74       endif
  75       if(.not.createvar(10,'z',1,LWORK,lDWORK)) return
  76
  77       EPS = dlamch('eps')
  78       ANORM = zlange( '1', M, N, zstk(lA), M, zstk(lDWORK) )
  79 c
  80       if(M.eq.N) then
  81 c
  82 c     M = N
  83 c
  84          call ZLACPY( 'F', N, N, zstk(lA), N, zstk(lAF), N )
  85 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
  86          call ZLACPY( 'F', N, NRHS, zstk(lB), N, zstk(lXB), N )
  87 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
  88          call ZGETRF( N, N, zstk(lAF), N, istk(lIPIV), INFO )
  89 c     SUBROUTINE ZGETRF( N, N, A, LDA, IPIV, INFO )
  90          if(INFO.eq.0) then
  91             call ZGECON( '1', N, zstk(lAF), N, ANORM, RCOND,
  92      $           zstk(lDWORK), stk(lRWORK), INFO )
  93 c     SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK,
  94 c     $                        RWORK, INFO )
  95             if(RCOND.gt.sqrt(EPS)) then
  96                call ZGETRS( 'N', N, NRHS, zstk(lAF), N, istk(lIPIV),
  97      $              zstk(lXB), N, INFO )
  98 c     SUBROUTINE ZGETRS( TRANS, N, NRHS, A, LDA, IPIV,
  99 c     B, LDB, INFO )
 100                call ZLACPY( 'F', N, NRHS, zstk(lXB), N, zstk(lX), N )
 101 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
 102                lhsvar(1) = 4
 103                return
 104             else
 105 c     .        ill conditioned problem
 106                call writebufzgesv3(buf,RCOND)
 107                call msgs(5,1)
 108             endif
 109          endif
 110       endif
 111 c
 112 c     M.ne.N or A  singular
 113 c
 114       RCOND = sqrt(EPS)
 115       call ZLACPY( 'F', M, NRHS, zstk(lB), M, zstk(lXB), max(M,N) )
 116 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
 117       do 10 i = 1, N
 118           istk(lJPVT+i-1) = 0
 119  10   continue
 120       call ZGELSY1( M, N, NRHS, zstk(lA), M, zstk(lXB), max(M,N),
 121      $     istk(lJPVT), RCOND, istk(lRANK), zstk(lDWORK),
 122      $     LWORK, stk(lRWORK), INFO )
 123 c     SUBROUTINE ZGELSY1( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
 124 c     $                    RANK, WORK, LWORK, RWORK, 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 ZLACPY( 'F', N, NRHS, zstk(lXB), max(M,N), zstk(lX), N )
 132 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
 133
 134       lhsvar(1)=4
 135       return
 136 c
 137       end
 138