sci2c arduino updated

1 files changed, 477 insertions, 0 deletions

diff --git a/2.3-1/src/fortran/lapack/zgges.f b/2.3-1/src/fortran/lapack/zgges.f
new file mode 100644
index 00000000..c1499003
--- /dev/null
+++ b/2.3-1/src/fortran/lapack/zgges.f
@@ -0,0 +1,477 @@
+      SUBROUTINE ZGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
+     $                  SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
+     $                  LWORK, RWORK, BWORK, INFO )
+*
+*  -- LAPACK driver routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVSL, JOBVSR, SORT
+      INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            BWORK( * )
+      DOUBLE PRECISION   RWORK( * )
+      COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
+     $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
+     $                   WORK( * )
+*     ..
+*     .. Function Arguments ..
+      LOGICAL            SELCTG
+      EXTERNAL           SELCTG
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZGGES computes for a pair of N-by-N complex nonsymmetric matrices
+*  (A,B), the generalized eigenvalues, the generalized complex Schur
+*  form (S, T), and optionally left and/or right Schur vectors (VSL
+*  and VSR). This gives the generalized Schur factorization
+*
+*          (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
+*
+*  where (VSR)**H is the conjugate-transpose of VSR.
+*
+*  Optionally, it also orders the eigenvalues so that a selected cluster
+*  of eigenvalues appears in the leading diagonal blocks of the upper
+*  triangular matrix S and the upper triangular matrix T. The leading
+*  columns of VSL and VSR then form an unitary basis for the
+*  corresponding left and right eigenspaces (deflating subspaces).
+*
+*  (If only the generalized eigenvalues are needed, use the driver
+*  ZGGEV instead, which is faster.)
+*
+*  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
+*  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
+*  usually represented as the pair (alpha,beta), as there is a
+*  reasonable interpretation for beta=0, and even for both being zero.
+*
+*  A pair of matrices (S,T) is in generalized complex Schur form if S
+*  and T are upper triangular and, in addition, the diagonal elements
+*  of T are non-negative real numbers.
+*
+*  Arguments
+*  =========
+*
+*  JOBVSL  (input) CHARACTER*1
+*          = 'N':  do not compute the left Schur vectors;
+*          = 'V':  compute the left Schur vectors.
+*
+*  JOBVSR  (input) CHARACTER*1
+*          = 'N':  do not compute the right Schur vectors;
+*          = 'V':  compute the right Schur vectors.
+*
+*  SORT    (input) CHARACTER*1
+*          Specifies whether or not to order the eigenvalues on the
+*          diagonal of the generalized Schur form.
+*          = 'N':  Eigenvalues are not ordered;
+*          = 'S':  Eigenvalues are ordered (see SELCTG).
+*
+*  SELCTG  (external procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
+*          SELCTG must be declared EXTERNAL in the calling subroutine.
+*          If SORT = 'N', SELCTG is not referenced.
+*          If SORT = 'S', SELCTG is used to select eigenvalues to sort
+*          to the top left of the Schur form.
+*          An eigenvalue ALPHA(j)/BETA(j) is selected if
+*          SELCTG(ALPHA(j),BETA(j)) is true.
+*
+*          Note that a selected complex eigenvalue may no longer satisfy
+*          SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
+*          ordering may change the value of complex eigenvalues
+*          (especially if the eigenvalue is ill-conditioned), in this
+*          case INFO is set to N+2 (See INFO below).
+*
+*  N       (input) INTEGER
+*          The order of the matrices A, B, VSL, and VSR.  N >= 0.
+*
+*  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
+*          On entry, the first of the pair of matrices.
+*          On exit, A has been overwritten by its generalized Schur
+*          form S.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of A.  LDA >= max(1,N).
+*
+*  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
+*          On entry, the second of the pair of matrices.
+*          On exit, B has been overwritten by its generalized Schur
+*          form T.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of B.  LDB >= max(1,N).
+*
+*  SDIM    (output) INTEGER
+*          If SORT = 'N', SDIM = 0.
+*          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
+*          for which SELCTG is true.
+*
+*  ALPHA   (output) COMPLEX*16 array, dimension (N)
+*  BETA    (output) COMPLEX*16 array, dimension (N)
+*          On exit,  ALPHA(j)/BETA(j), j=1,...,N, will be the
+*          generalized eigenvalues.  ALPHA(j), j=1,...,N  and  BETA(j),
+*          j=1,...,N  are the diagonals of the complex Schur form (A,B)
+*          output by ZGGES. The  BETA(j) will be non-negative real.
+*
+*          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
+*          underflow, and BETA(j) may even be zero.  Thus, the user
+*          should avoid naively computing the ratio alpha/beta.
+*          However, ALPHA will be always less than and usually
+*          comparable with norm(A) in magnitude, and BETA always less
+*          than and usually comparable with norm(B).
+*
+*  VSL     (output) COMPLEX*16 array, dimension (LDVSL,N)
+*          If JOBVSL = 'V', VSL will contain the left Schur vectors.
+*          Not referenced if JOBVSL = 'N'.
+*
+*  LDVSL   (input) INTEGER
+*          The leading dimension of the matrix VSL. LDVSL >= 1, and
+*          if JOBVSL = 'V', LDVSL >= N.
+*
+*  VSR     (output) COMPLEX*16 array, dimension (LDVSR,N)
+*          If JOBVSR = 'V', VSR will contain the right Schur vectors.
+*          Not referenced if JOBVSR = 'N'.
+*
+*  LDVSR   (input) INTEGER
+*          The leading dimension of the matrix VSR. LDVSR >= 1, and
+*          if JOBVSR = 'V', LDVSR >= N.
+*
+*  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
+*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK.  LWORK >= max(1,2*N).
+*          For good performance, LWORK must generally be larger.
+*
+*          If LWORK = -1, then a workspace query is assumed; the routine
+*          only calculates the optimal size of the WORK array, returns
+*          this value as the first entry of the WORK array, and no error
+*          message related to LWORK is issued by XERBLA.
+*
+*  RWORK   (workspace) DOUBLE PRECISION array, dimension (8*N)
+*
+*  BWORK   (workspace) LOGICAL array, dimension (N)
+*          Not referenced if SORT = 'N'.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value.
+*          =1,...,N:
+*                The QZ iteration failed.  (A,B) are not in Schur
+*                form, but ALPHA(j) and BETA(j) should be correct for
+*                j=INFO+1,...,N.
+*          > N:  =N+1: other than QZ iteration failed in ZHGEQZ
+*                =N+2: after reordering, roundoff changed values of
+*                      some complex eigenvalues so that leading
+*                      eigenvalues in the Generalized Schur form no
+*                      longer satisfy SELCTG=.TRUE.  This could also
+*                      be caused due to scaling.
+*                =N+3: reordering falied in ZTGSEN.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+      COMPLEX*16         CZERO, CONE
+      PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
+     $                   CONE = ( 1.0D0, 0.0D0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
+     $                   LQUERY, WANTST
+      INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
+     $                   ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, LWKMIN,
+     $                   LWKOPT
+      DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
+     $                   PVSR, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      INTEGER            IDUM( 1 )
+      DOUBLE PRECISION   DIF( 2 )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
+     $                   ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGSEN, ZUNGQR,
+     $                   ZUNMQR
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      DOUBLE PRECISION   DLAMCH, ZLANGE
+      EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+*     Decode the input arguments
+*
+      IF( LSAME( JOBVSL, 'N' ) ) THEN
+         IJOBVL = 1
+         ILVSL = .FALSE.
+      ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
+         IJOBVL = 2
+         ILVSL = .TRUE.
+      ELSE
+         IJOBVL = -1
+         ILVSL = .FALSE.
+      END IF
+*
+      IF( LSAME( JOBVSR, 'N' ) ) THEN
+         IJOBVR = 1
+         ILVSR = .FALSE.
+      ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
+         IJOBVR = 2
+         ILVSR = .TRUE.
+      ELSE
+         IJOBVR = -1
+         ILVSR = .FALSE.
+      END IF
+*
+      WANTST = LSAME( SORT, 'S' )
+*
+*     Test the input arguments
+*
+      INFO = 0
+      LQUERY = ( LWORK.EQ.-1 )
+      IF( IJOBVL.LE.0 ) THEN
+         INFO = -1
+      ELSE IF( IJOBVR.LE.0 ) THEN
+         INFO = -2
+      ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -7
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
+         INFO = -14
+      ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
+         INFO = -16
+      END IF
+*
+*     Compute workspace
+*      (Note: Comments in the code beginning "Workspace:" describe the
+*       minimal amount of workspace needed at that point in the code,
+*       as well as the preferred amount for good performance.
+*       NB refers to the optimal block size for the immediately
+*       following subroutine, as returned by ILAENV.)
+*
+      IF( INFO.EQ.0 ) THEN
+         LWKMIN = MAX( 1, 2*N )
+         LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
+         LWKOPT = MAX( LWKOPT, N +
+     $                 N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, -1 ) )
+         IF( ILVSL ) THEN
+            LWKOPT = MAX( LWKOPT, N +
+     $                    N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
+         END IF
+         WORK( 1 ) = LWKOPT
+*
+         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
+     $      INFO = -18
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZGGES ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 ) THEN
+         SDIM = 0
+         RETURN
+      END IF
+*
+*     Get machine constants
+*
+      EPS = DLAMCH( 'P' )
+      SMLNUM = DLAMCH( 'S' )
+      BIGNUM = ONE / SMLNUM
+      CALL DLABAD( SMLNUM, BIGNUM )
+      SMLNUM = SQRT( SMLNUM ) / EPS
+      BIGNUM = ONE / SMLNUM
+*
+*     Scale A if max element outside range [SMLNUM,BIGNUM]
+*
+      ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
+      ILASCL = .FALSE.
+      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
+         ANRMTO = SMLNUM
+         ILASCL = .TRUE.
+      ELSE IF( ANRM.GT.BIGNUM ) THEN
+         ANRMTO = BIGNUM
+         ILASCL = .TRUE.
+      END IF
+*
+      IF( ILASCL )
+     $   CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
+*
+*     Scale B if max element outside range [SMLNUM,BIGNUM]
+*
+      BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
+      ILBSCL = .FALSE.
+      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
+         BNRMTO = SMLNUM
+         ILBSCL = .TRUE.
+      ELSE IF( BNRM.GT.BIGNUM ) THEN
+         BNRMTO = BIGNUM
+         ILBSCL = .TRUE.
+      END IF
+*
+      IF( ILBSCL )
+     $   CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
+*
+*     Permute the matrix to make it more nearly triangular
+*     (Real Workspace: need 6*N)
+*
+      ILEFT = 1
+      IRIGHT = N + 1
+      IRWRK = IRIGHT + N
+      CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
+     $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
+*
+*     Reduce B to triangular form (QR decomposition of B)
+*     (Complex Workspace: need N, prefer N*NB)
+*
+      IROWS = IHI + 1 - ILO
+      ICOLS = N + 1 - ILO
+      ITAU = 1
+      IWRK = ITAU + IROWS
+      CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
+     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
+*
+*     Apply the orthogonal transformation to matrix A
+*     (Complex Workspace: need N, prefer N*NB)
+*
+      CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
+     $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
+     $             LWORK+1-IWRK, IERR )
+*
+*     Initialize VSL
+*     (Complex Workspace: need N, prefer N*NB)
+*
+      IF( ILVSL ) THEN
+         CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
+         IF( IROWS.GT.1 ) THEN
+            CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                   VSL( ILO+1, ILO ), LDVSL )
+         END IF
+         CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
+     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
+      END IF
+*
+*     Initialize VSR
+*
+      IF( ILVSR )
+     $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
+*
+*     Reduce to generalized Hessenberg form
+*     (Workspace: none needed)
+*
+      CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
+     $             LDVSL, VSR, LDVSR, IERR )
+*
+      SDIM = 0
+*
+*     Perform QZ algorithm, computing Schur vectors if desired
+*     (Complex Workspace: need N)
+*     (Real Workspace: need N)
+*
+      IWRK = ITAU
+      CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
+     $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ),
+     $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
+      IF( IERR.NE.0 ) THEN
+         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
+            INFO = IERR
+         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
+            INFO = IERR - N
+         ELSE
+            INFO = N + 1
+         END IF
+         GO TO 30
+      END IF
+*
+*     Sort eigenvalues ALPHA/BETA if desired
+*     (Workspace: none needed)
+*
+      IF( WANTST ) THEN
+*
+*        Undo scaling on eigenvalues before selecting
+*
+         IF( ILASCL )
+     $      CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, 1, ALPHA, N, IERR )
+         IF( ILBSCL )
+     $      CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, 1, BETA, N, IERR )
+*
+*        Select eigenvalues
+*
+         DO 10 I = 1, N
+            BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) )
+   10    CONTINUE
+*
+         CALL ZTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHA,
+     $                BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL, PVSR,
+     $                DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1, IERR )
+         IF( IERR.EQ.1 )
+     $      INFO = N + 3
+*
+      END IF
+*
+*     Apply back-permutation to VSL and VSR
+*     (Workspace: none needed)
+*
+      IF( ILVSL )
+     $   CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
+     $                RWORK( IRIGHT ), N, VSL, LDVSL, IERR )
+      IF( ILVSR )
+     $   CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
+     $                RWORK( IRIGHT ), N, VSR, LDVSR, IERR )
+*
+*     Undo scaling
+*
+      IF( ILASCL ) THEN
+         CALL ZLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
+         CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
+      END IF
+*
+      IF( ILBSCL ) THEN
+         CALL ZLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
+         CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
+      END IF
+*
+      IF( WANTST ) THEN
+*
+*        Check if reordering is correct
+*
+         LASTSL = .TRUE.
+         SDIM = 0
+         DO 20 I = 1, N
+            CURSL = SELCTG( ALPHA( I ), BETA( I ) )
+            IF( CURSL )
+     $         SDIM = SDIM + 1
+            IF( CURSL .AND. .NOT.LASTSL )
+     $         INFO = N + 2
+            LASTSL = CURSL
+   20    CONTINUE
+*
+      END IF
+*
+   30 CONTINUE
+*
+      WORK( 1 ) = LWKOPT
+*
+      RETURN
+*
+*     End of ZGGES
+*
+      END