MINICERN/mathlib/gen/e/decomp.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:27  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE DECOMP(M,N,QR,ALPHA,IPIVOT,ERR,YA,SUM)
  11       DIMENSION QR(M,N),ALPHA(N),IPIVOT(N),YA(N),SUM(N)
  12       LOGICAL ERR
  13 C
  14 C        THIS ROUTINE REDUCES THE MATRIX GIVEN IN THE ARRAY
  15 C        QR(M,N), WHERE M \ N, TO UPPER RIGHT TRIANGULAR FORM
  16 C        BY MEANS OF N ELEMENTARY ORTHOGONAL TRANSFORMATIONS
  17 C        ( I - BETA UT , WHERE T IS THE TRANSPOSE OF U).
  18 C        THE DIAGONAL ELEMENTS OF THE REDUCED MATRIX ARE STORED
  19 C        IN THE ARRAY ALPHA(N), THE OFF-DIAGONAL ELEMENTS IN THE
  20 C        UPPER RIGHT TRIANGULAR PART OF QR.  THE NON-ZERO
  21 C        COMPONENTS OF THE VECTORS U ARE STORED ON AND BELOW
  22 C        THE LEADING DIAGONAL OF QR.  PIVOTING IS DONE BY
  23 C        CHOOSING AT EACH STEP, THE COLUMN WITH THE LARGEST SUM
  24 C        OF SQUARES TO BE REDUCED NEXT.  THESE INTERCHANGES ARE
  25 C        RECORDED IN THE ARRAY IPIVOT(N).  IF AT ANY STAGE, THE
  26 C        SUM OF SQUARES OF THE COLUMN TO BE REDUCED IS EXACTLY
  27 C        EQUAL TO ZERO, THEN THE LOGICAL VARIABLE ERR IS SET TO
  28 C        .FALSE. AND CONTROL IS RETURNED TO THE MAIN PROGRAM.
  29 C
  30       DO 5  J=1,N
  31 C     ***  J TH. COLUMN SUM  ***
  32       SUM(J)=PROD1(QR(1,J),QR(1,J),1,M)
  33  5    IPIVOT(J)=J
  34 C
  35       DO 40  K=1,N
  36 C     ***  K TH. HOUSEHOLDER TRANSFORMATION  ***
  37       SIGMA=SUM(K)
  38       JBAR=K
  39       L=K+1
  40       IF(L.GT.N) GO TO 11
  41       DO 10  J=L,N
  42       IF(SIGMA .GE. SUM(J)) GO TO 10
  43       SIGMA=SUM(J)
  44       JBAR=J
  45  10   CONTINUE
  46  11   CONTINUE
  47       IF(JBAR .EQ. K) GO TO 20
  48 C     ***  COLUMN INTERCHANGE  ***
  49       I=IPIVOT(K)
  50       IPIVOT(K)=IPIVOT(JBAR)
  51       IPIVOT(JBAR)=I
  52       SUM(JBAR)=SUM(K)
  53       SUM(K)=SIGMA
  54       DO 15  I=1,M
  55       SIGMA=QR(I,K)
  56       QR(I,K)=QR(I,JBAR)
  57  15   QR(I,JBAR)=SIGMA
  58  20   CONTINUE
  59       SIGMA=PROD1(QR(1,K),QR(1,K),K,M)
  60       IF(SIGMA .NE. 0.0) GO TO 22
  61       ERR=.FALSE.
  62       RETURN
  63  22   QRKK=QR(K,K)
  64       ALPHA(K)=SQRT(SIGMA)
  65       IF(QRKK .GE. 0.0)  ALPHA(K)=-ALPHA(K)
  66       ALPHAK=ALPHA(K)
  67       BETA=1.0/(SIGMA-QRKK*ALPHAK)
  68       QR(K,K)=QRKK-ALPHAK
  69       IF (L.GT.N) GO TO 40
  70       DO 25  J=L,N
  71  25   YA(J)=BETA*PROD1(QR(1,K),QR(1,J),K,M)
  72       DO 35  J=L,N
  73       DO 30  I=K,M
  74  30   QR(I,J)=QR(I,J)-QR(I,K)*YA(J)
  75  35   SUM(J)=SUM(J)-QR(K,J)**2
  76  40   CONTINUE
  77       RETURN
  78       END