669c72e94b6ee320226e95e13c5b0d378919b818
[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$
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