*
* $Id: decomp.F,v 1.1.1.1 1996/04/01 15:02:27 mclareni Exp $
*
* $Log: decomp.F,v $
* Revision 1.1.1.1  1996/04/01 15:02:27  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE DECOMP(M,N,QR,ALPHA,IPIVOT,ERR,YA,SUM)
      DIMENSION QR(M,N),ALPHA(N),IPIVOT(N),YA(N),SUM(N)
      LOGICAL ERR
C
C        THIS ROUTINE REDUCES THE MATRIX GIVEN IN THE ARRAY
C        QR(M,N), WHERE M \ N, TO UPPER RIGHT TRIANGULAR FORM
C        BY MEANS OF N ELEMENTARY ORTHOGONAL TRANSFORMATIONS
C        ( I - BETA UT , WHERE T IS THE TRANSPOSE OF U).
C        THE DIAGONAL ELEMENTS OF THE REDUCED MATRIX ARE STORED
C        IN THE ARRAY ALPHA(N), THE OFF-DIAGONAL ELEMENTS IN THE
C        UPPER RIGHT TRIANGULAR PART OF QR.  THE NON-ZERO
C        COMPONENTS OF THE VECTORS U ARE STORED ON AND BELOW
C        THE LEADING DIAGONAL OF QR.  PIVOTING IS DONE BY
C        CHOOSING AT EACH STEP, THE COLUMN WITH THE LARGEST SUM
C        OF SQUARES TO BE REDUCED NEXT.  THESE INTERCHANGES ARE
C        RECORDED IN THE ARRAY IPIVOT(N).  IF AT ANY STAGE, THE
C        SUM OF SQUARES OF THE COLUMN TO BE REDUCED IS EXACTLY
C        EQUAL TO ZERO, THEN THE LOGICAL VARIABLE ERR IS SET TO
C        .FALSE. AND CONTROL IS RETURNED TO THE MAIN PROGRAM.
C
      DO 5  J=1,N
C     ***  J TH. COLUMN SUM  ***
      SUM(J)=PROD1(QR(1,J),QR(1,J),1,M)
 5    IPIVOT(J)=J
C
      DO 40  K=1,N
C     ***  K TH. HOUSEHOLDER TRANSFORMATION  ***
      SIGMA=SUM(K)
      JBAR=K
      L=K+1
      IF(L.GT.N) GO TO 11
      DO 10  J=L,N
      IF(SIGMA .GE. SUM(J)) GO TO 10
      SIGMA=SUM(J)
      JBAR=J
 10   CONTINUE
 11   CONTINUE
      IF(JBAR .EQ. K) GO TO 20
C     ***  COLUMN INTERCHANGE  ***
      I=IPIVOT(K)
      IPIVOT(K)=IPIVOT(JBAR)
      IPIVOT(JBAR)=I
      SUM(JBAR)=SUM(K)
      SUM(K)=SIGMA
      DO 15  I=1,M
      SIGMA=QR(I,K)
      QR(I,K)=QR(I,JBAR)
 15   QR(I,JBAR)=SIGMA
 20   CONTINUE
      SIGMA=PROD1(QR(1,K),QR(1,K),K,M)
      IF(SIGMA .NE. 0.0) GO TO 22
      ERR=.FALSE.
      RETURN
 22   QRKK=QR(K,K)
      ALPHA(K)=SQRT(SIGMA)
      IF(QRKK .GE. 0.0)  ALPHA(K)=-ALPHA(K)
      ALPHAK=ALPHA(K)
      BETA=1.0/(SIGMA-QRKK*ALPHAK)
      QR(K,K)=QRKK-ALPHAK
      IF (L.GT.N) GO TO 40
      DO 25  J=L,N
 25   YA(J)=BETA*PROD1(QR(1,K),QR(1,J),K,M)
      DO 35  J=L,N
      DO 30  I=K,M
 30   QR(I,J)=QR(I,J)-QR(I,K)*YA(J)
 35   SUM(J)=SUM(J)-QR(K,J)**2
 40   CONTINUE
      RETURN
      END