* Bug #9196 fixed - The threshold level for conditioning in backslash was too
[scilab.git] / 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       double precision RCOND_thresh
22       integer vfinite
23       external dlamch, zlange, vfinite
24       intrinsic sqrt
25 c     
26       minrhs=2
27       maxrhs=2
28       minlhs=1
29       maxlhs=1
30
31 c     
32       if(.not.checkrhs(fname,minrhs,maxrhs)) return
33       if(.not.checklhs(fname,minlhs,maxlhs)) return
34       
35       if(.not.getrhsvar(1,'z', MA, NA, lA)) return
36       if(.not.getrhsvar(2,'z', 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,'z', 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*2, zstk(lA)).eq.0) then
54          call error(229)
55          return
56       endif
57       if (vfinite(MB*NRHS*2, zstk(lB)).eq.0) then
58          call error(229)
59          return
60       endif
61
62       if(.not.createvar(3,'z', M, N, lAF)) return
63       if(.not.createvar(4,'z', N, NRHS, lX)) return
64       if(.not.createvar(5,'z', max(M,N), NRHS, lXB)) return
65       if(.not.createvar(6,'i', 1, 1, lRANK)) return
66       if(.not.createvar(7,'i', 1, N, lIPIV)) return
67       if(.not.createvar(8,'i', 1, N, lJPVT)) return
68       if(.not.createvar(9,'d',1,2*N,lRWORK)) return       
69       LWORKMIN = max( 2*N, min(M,N)+max( 2*min(M,N), N+1,
70      $     min(M,N)+NRHS ) )
71       LWORK=maxvol(10,'z')
72       if(LWORK.le.LWORKMIN) then
73          err=2*(LWORK-LWORKMIN)
74          call error(17)
75          return
76       endif
77       if(.not.createvar(10,'z',1,LWORK,lDWORK)) return
78       
79       EPS = dlamch('eps')
80       RCOND_thresh=EPS*10
81       ANORM = zlange( '1', M, N, zstk(lA), M, zstk(lDWORK) )
82 c     
83       if(M.eq.N) then
84 c     
85 c     M = N
86 c     
87          call ZLACPY( 'F', N, N, zstk(lA), N, zstk(lAF), N )
88 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
89          call ZLACPY( 'F', N, NRHS, zstk(lB), N, zstk(lXB), N )
90 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
91          call ZGETRF( N, N, zstk(lAF), N, istk(lIPIV), INFO )         
92 c     SUBROUTINE ZGETRF( N, N, A, LDA, IPIV, INFO )
93          if(INFO.eq.0) then
94             call ZGECON( '1', N, zstk(lAF), N, ANORM, RCOND,
95      $           zstk(lDWORK), stk(lRWORK), INFO )
96 c     SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK,
97 c     $                        RWORK, INFO )
98             if(RCOND.gt.RCOND_thresh) then
99                call ZGETRS( 'N', N, NRHS, zstk(lAF), N, istk(lIPIV),
100      $              zstk(lXB), N, INFO ) 
101 c     SUBROUTINE ZGETRS( TRANS, N, NRHS, A, LDA, IPIV,
102 c     B, LDB, INFO )
103                call ZLACPY( 'F', N, NRHS, zstk(lXB), N, zstk(lX), N )
104 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
105                lhsvar(1) = 4
106                return
107             else
108 c     .        ill conditioned problem
109                call writebufzgesv3(buf,RCOND)
110                call msgs(5,1)
111             endif
112          endif
113       endif
114 c     
115 c     M.ne.N or A  singular
116 c     
117       RCOND = RCOND_thresh
118       call ZLACPY( 'F', M, NRHS, zstk(lB), M, zstk(lXB), max(M,N) )
119 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
120       do 10 i = 1, N
121           istk(lJPVT+i-1) = 0
122  10   continue
123       call ZGELSY1( M, N, NRHS, zstk(lA), M, zstk(lXB), max(M,N),
124      $     istk(lJPVT), RCOND, istk(lRANK), zstk(lDWORK),
125      $     LWORK, stk(lRWORK), INFO )
126 c     SUBROUTINE ZGELSY1( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
127 c     $                    RANK, WORK, LWORK, RWORK, INFO )
128       if(info.ne.0) then
129          return
130       endif
131       if( M.ne.N .and. istk(lRANK).lt.min(M,N) )then
132          call msgs(9,istk(lRANK))
133       endif
134       call ZLACPY( 'F', N, NRHS, zstk(lXB), max(M,N), zstk(lX), N )
135 c     SUBROUTINE ZLACPY( UPLO, M, N, A, LDA, B, LDB )
136
137       lhsvar(1)=4
138       return
139 c     
140       end
141