[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / decomp.F

*
* $Id$
*
* $Log$
* 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
Commit	Line	Data
fe4da5cc	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