* Bug #9196 fixed - The threshold level for conditioning in backslash was too
[scilab.git] / 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