6ceb4ce65a901ea45a7ba306b175f1eabbddf788
[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       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