C HAPACK_SYMPQR.F
C   Gateway function for the symplectic QR factorization of a
C   general matrix.
C
C Matlab call:  
C              [(Q1,Q2,)R] = hapack_sympqr(task,A)
C
C   task =  1 :        [R] = hapack_sympqr(1,A)
C                [Q1,Q2,R] = hapack_sympqr(1,A)
C
C   task =  2 :        [R] = hapack_sympqr(2,A)
C                [Q1,Q2,R] = hapack_sympqr(2,A)
C
C Purpose:
C   To compute a symplectic QR factorization of a general 2*M-by-N
C   matrix A.
C
C   task =  1: Compute the R factor R = [R1; R2], where R1 is an M-by-N
C              upper trapezoidal matrix and R2 is an M-by-N strictly
C              upper trapezoidal matrix. Optionally, M-by-M matrices
C              Q1 and Q2 are computed so that Q = [Q1, Q2; -Q2, Q1]
C              satisfies Q'*Q=I and A = Q*R.
C
C   task =  2: Compute R and "economy size" matrices Q1 and Q2 of, i.e.,
C              only their first min(M,N) columns are returned.
C
C Input parameters: 
C   task  - integer option to determine the computation to perform as
C           described above.
C   A     - real 2*M-by-N matrix.
C
C Output parameters:
C   Q1,Q2 - real M-by-K matrices so that Q = [Q1, Q2; -Q2, Q1]
C           satisfies Q'*Q=I and A = Q*R. If task = 1 then K = M,
C           and if task = 2 then K = MIN(M,N).
C   R     - the 2*M-by-N upper trapezoidal/strictly upper trapezoidal
C           factor.
C
C Contributor:
C   D. Kressner, UCL CESAME, Louvain-la-Neuve, March 2003.
C
C **********************************************************************
C
      SUBROUTINE MEXFUNCTION( NLHS, PLHS, NRHS, PRHS )
      IMPLICIT NONE
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO
      PARAMETER         ( ZERO = 0.0D0 )
C
C     .. Mex-file interface parameters ..
      INTEGER           PLHS(*), PRHS(*)
      INTEGER           NLHS, NRHS
C
C     .. Mex-file integer functions ..
      INTEGER           mxCreateFull, mxGetPr
      INTEGER           mxGetM, mxGetN, mxIsComplex, mxIsNumeric
C
C     .. Scalar parameters used by HAPACK subroutines ..
      INTEGER           INFO, LDA, LDWORK, LDQ1, LDQ2, M, N
C
C     .. Allocatable arrays ..
C     !Fortran 90/95 (Fixed dimensions should be used with Fortran 77.)
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:), CS(:), DWORK(:), Q1(:,:),
     $                                 Q2(:,:), TAU(:)
C
C     .. Local variables and constant dimension arrays ..
      CHARACTER*120     TEXT
      LOGICAL           COMPQ
      INTEGER           IP, K, MM, MN, NB, NBMIN, TASK
      DOUBLE PRECISION  TEMP
C
C     .. External Functions ..
      INTEGER           ILAHAP
      EXTERNAL          ILAHAP
C
C     .. External subroutines ..
      EXTERNAL          DGESQB, DLACPY, DLASET, DOSGSB
C
C     ..Intrinsic functions..
      INTRINSIC         MAX, MIN
C
C     Check for proper number of arguments.
C
      IF ( NRHS.NE.2 ) THEN
         CALL mexErrMsgTxt
     $        ( 'HAPACK_SYMPQR requires 2 input arguments' )
      ELSE IF ( NLHS.GT.3 ) THEN
         CALL mexErrMsgTxt
     $        ( 'HAPACK_SYMPQR requires at most 3 output arguments' )
      END IF
      IF ( NLHS.GT.1 ) THEN
         COMPQ = .TRUE.
      ELSE
         COMPQ = .FALSE.
      END IF
C
C     Check dimensions of input parameters and read/set scalar
C     parameters.
C
      IF ( mxGetM( PRHS(1) ).NE.1 .OR. mxGetN( PRHS(1) ).NE.1 )
     $   CALL mexErrMsgTxt( 'TASK must be a scalar' )
      IF ( mxIsNumeric( PRHS(1) ).EQ.0 .OR. 
     $     mxIsComplex( PRHS(1) ).EQ.1 )
     $   CALL mexErrMsgTxt( 'TASK must be an integer scalar' )
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(1) ), TEMP, 1 )
      TASK = TEMP
C
      IF ( TASK.LT.1 .OR. TASK.GT.2 )
     $   CALL mexErrMsgTxt
     $        ( 'The only admissible values of TASK are 1 and 2')
C
C     Check A.
C
      MM = mxGetM( PRHS(2) )
      N = mxGetN( PRHS(2) )
      IF ( MM.LT.0 .OR. N.LT.0 )
     $   CALL mexErrMsgTxt( 'A must be a matrix' )
      M = MM / 2
      IF ( M+M.NE.MM )
     $   CALL mexErrMsgTxt( 'A must have an even number of rows' )
      IF ( mxIsNumeric( PRHS(2) ).EQ.0 .OR. 
     $     mxIsComplex( PRHS(2) ).EQ.1 )
     $   CALL mexErrMsgTxt( 'A must be a real matrix' )
C
      MN = MIN( M, N )
      IF ( TASK.EQ.1 ) THEN
         K = M
      ELSE
         K = MN
      END IF
C
C     Determine leading dimensions.
C
      LDA = MAX( 1, MM )
      IF ( COMPQ ) THEN
         LDQ1 = MAX( 1, M )
         LDQ2 = MAX( 1, M )
      END IF
C
C     Determine length of workspace to enable blocked algorithms.
C
      LDWORK = N
      NB = ILAHAP( 1, 'DGESQB', '', M, N, -1 )
      NBMIN = MAX( 2, ILAHAP( 2, 'DGESQB',  ' ', M, N, -1 ) )
      IF ( NB.GE.NBMIN )
     $   LDWORK = MAX( LDWORK, 9*N*NB + 15*NB*NB )
      IF ( COMPQ ) THEN
         LDWORK = MAX( LDWORK, M+N )
         NB = ILAHAP( 1, 'DOSGSB', 'NN', M, N, MN )
         NBMIN = MAX( 2, ILAHAP( 2, 'DOSGSB', 'NN', M, N, MN ) )
         IF ( NB.GE.NBMIN )
     $      LDWORK = MAX( LDWORK, 8*N*NB + 15*NB*NB )
      END IF
C
C     Allocate variable dimension local arrays.
C     !Fortran 90/95
C
      ALLOCATE ( A(LDA,N), CS(2*MN), DWORK(LDWORK), TAU(MN) )
      IF ( COMPQ )
     $   ALLOCATE( Q1(LDQ1,K), Q2(LDQ2,K) )
C
C     Copy inputs from MATLAB workspace to locally allocated arrays.
C
      CALL mxCopyPtrToReal8( mxGetPr( PRHS(2) ), A, MM*N )
C
C     Do the actual computations.
C
      CALL DGESQB( M, N, A, LDA, A(M+1,1), LDA, CS, TAU, DWORK, LDWORK,
     $             INFO )
      IF ( INFO.EQ.0 ) THEN
         IF ( COMPQ ) THEN
            IF ( M.GT.1 )
     $         CALL DLACPY( 'Lower', M-1, MN, A(2,1), LDA, Q1(2,1),
     $                      LDQ1 )
            CALL DLACPY( 'Lower', M, MN, A(M+1,1), LDA, Q2, LDQ2 )
            CALL DOSGSB( 'No Transpose', 'No Transpose', M, K, MN, Q1,
     $                   LDQ1, Q2, LDQ2, CS, TAU, DWORK, LDWORK, INFO )
         END IF
         CALL DLASET( 'Lower', M-1, MIN(M,N), ZERO, ZERO, A(2,1), LDA )
         CALL DLASET( 'Lower', M, MIN(M,N), ZERO, ZERO, A(M+1,1), LDA )
      END IF
C
C     Copy output to MATLAB workspace.
C
      IP = 1
      IF ( INFO.EQ.0 ) THEN
         IF ( COMPQ ) THEN
            PLHS(IP) = mxCreateFull( M, K, 0 )
            CALL mxCopyReal8ToPtr( Q1, mxGetPr( PLHS(IP) ), M*K )
            IP = IP + 1
            PLHS(IP) = mxCreateFull( M, K, 0 )
            CALL mxCopyReal8ToPtr( Q2, mxGetPr( PLHS(IP) ), M*K )
            IP = IP + 1
         END IF
         IF ( NLHS.GE.IP ) THEN
            PLHS(IP) = mxCreateFull( MM, N, 0 )
            CALL mxCopyReal8ToPtr( A, mxGetPr( PLHS(IP) ), MM*N )
            IP = IP + 1
         END IF
      END IF
C
C     Deallocate local arrays.
C     !Fortran 90/95
C
      DEALLOCATE( A, CS, DWORK, TAU )
      IF ( COMPQ )
     $   DEALLOCATE( Q1, Q2 )
C
C     Error and warning handling.
C
      IF ( INFO.NE.0 ) THEN
         IF ( COMPQ ) THEN
            WRITE( TEXT, '('' INFO = '',I4,
     $             '' ON EXIT FROM DGESQB/DOSGSB'')' )  INFO
         ELSE
            WRITE( TEXT, '('' INFO = '',I4,'' ON EXIT FROM DGESQB'')' )
     $             INFO
         END IF
         CALL mexErrMsgTxt( TEXT )
      END IF
C
      RETURN
C *** Last line of HAPACK_SYMPQR ***
      END