Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dmchb.doc
There are 2 other files named dmchb.doc in the archive. Click here to see a list.
SUBROUTINE DMCHB

PURPOSE
   FOR A GIVEN POSITIVE-DEFINITE M BY M MATRIX A WITH SYMMETRIC
   BAND STRUCTURE AND - IF NECESSARY - A GIVEN GENERAL M BY N
   MATRIX R, THE FOLLOWING CALCULATIONS (DEPENDENT ON THE
   VALUE OF THE DECISION PARAMETER IOP) ARE PERFORMED
   (1) MATRIX A IS FACTORIZED (IF IOP IS NOT NEGATIVE), THAT
       MEANS BAND MATRIX TU WITH UPPER CODIAGONALS ONLY IS
       GENERATED ON THE LOCATIONS OF A SUCH THAT
       TRANSPOSE(TU)*TU=A.
   (2) MATRIX R IS MULTIPLIED ON THE LEFT BY INVERSE(TU)
       AND/OR INVERSE(TRANSPOSE(TU)) AND THE RESULT IS STORED
       IN THE LOCATIONS OF R.
   THIS SUBROUTINE ESPECIALLY CAN BE USED TO SOLVE THE SYSTEM
   OF SIMULTANEOUS LINEAR EQUATIONS A*X=R WITH POSITIVE-
   DEFINITE COEFFICIENT MATRIX A OF SYMMETRIC BAND STRUCTURE.

USAGE
   CALL DMCHB (R,A,M,N,MUD,IOP,EPS,IER)

DESCRIPTION OF PARAMETERS
   R	  - INPUT IN CASES IOP=-3,-2,-1,1,2,3  DOUBLE PRECISION
		  M BY N RIGHT HAND SIDE MATRIX,
		  IN CASE IOP=0  IRRELEVANT.
	    OUTPUT IN CASES IOP=1,-1  INVERSE(A)*R,
		   IN CASES IOP=2,-2  INVERSE(TU)*R,
		   IN CASES IOP=3,-3  INVERSE(TRANSPOSE(TU))*R,
		   IN CASE  IOP=0     UNCHANGED.
   A	  - INPUT IN CASES IOP=0,1,2,3	DOUBLE PRECISION M BY M
		  POSITIVE-DEFINITE COEFFICIENT MATRIX OF
		  SYMMETRIC BAND STRUCTURE STORED IN
		  COMPRESSED FORM (SEE REMARKS),
		  IN CASES IOP=-1,-2,-3 DOUBLE PRECISION M BY M
		  BAND MATRIX TU WITH UPPER CODIAGONALS ONLY,
		  STORED IN COMPRESSED FORM (SEE REMARKS).
	    OUTPUT IN ALL CASES  BAND MATRIX TU WITH UPPER
		   CODIAGONALS ONLY, STORED IN COMPRESSED FORM
		   (THAT MEANS UNCHANGED IF IOP=-1,-2,-3).
   M	  - INPUT VALUE SPECIFYING THE NUMBER OF ROWS AND
	    COLUMNS OF A AND THE NUMBER OF ROWS OF R.
   N	  - INPUT VALUE SPECIFYING THE NUMBER OF COLUMNS OF R
	    (IRRELEVANT IN CASE IOP=0).
   MUD	  - INPUT VALUE SPECIFYING THE NUMBER OF UPPER
	    CODIAGONALS OF A.
   IOP	  - ONE OF THE VALUES -3,-2,-1,0,1,2,3 GIVEN AS INPUT
	    AND USED AS DECISION PARAMETER.
   EPS	  - SINGLE PRECISION INPUT VALUE USED AS RELATIVE
	    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANT DIGITS.
   IER	  - RESULTING ERROR PARAMETER CODED AS FOLLOWS
	     IER=0  - NO ERROR,
	     IER=-1 - NO RESULT BECAUSE OF WRONG INPUT
		      PARAMETERS M,MUD,IOP (SEE REMARKS),
		      OR BECAUSE OF A NONPOSITIVE RADICAND AT
		      SOME FACTORIZATION STEP,
		      OR BECAUSE OF A ZERO DIAGONAL ELEMENT
		      AT SOME DIVISION STEP.
	     IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
		      CANCE INDICATED AT FACTORIZATION STEP K+1
		      WHERE RADICAND WAS NO LONGER GREATER
		      THAN EPS*A(K+1,K+1).

REMARKS
   UPPER PART OF SYMMETRIC BAND MATRIX A CONSISTING OF MAIN
   DIAGONAL AND MUD UPPER CODIAGONALS (RESP. BAND MATRIX TU
   CONSISTING OF MAIN DIAGONAL AND MUD UPPER CODIAGONALS)
   IS ASSUMED TO BE STORED IN COMPRESSED FORM, I.E. ROWWISE
   IN TOTALLY NEEDED M+MUD*(2M-MUD-1)/2 SUCCESSIVE STORAGE
   LOCATIONS. ON RETURN UPPER BAND FACTOR TU (ON THE LOCATIONS
   OF A) IS STORED IN THE SAME WAY.
   RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
   IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN RESULT MATRIX
   INVERSE(A)*R OR INVERSE(TU)*R OR INVERSE(TRANSPOSE(TU))*R
   IS STORED COLUMNWISE TOO ON THE LOCATIONS OF R.
   INPUT PARAMETERS M, MUD, IOP SHOULD SATISFY THE FOLLOWING
   RESTRICTIONS     MUD NOT LESS THAN ZERO,
		    1+MUD NOT GREATER THAN M,
		    ABS(IOP) NOT GREATER THAN 3.
   NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
   RESTRICTIONS ARE NOT SATISFIED.
   THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
   PARAMETERS ARE SATISFIED, IF RADICANDS AT ALL FACTORIZATION
   STEPS ARE POSITIVE AND/OR IF ALL DIAGONAL ELEMENTS OF
   UPPER BAND FACTOR TU ARE NONZERO.

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
   NONE

METHOD
   FACTORIZATION IS DONE USING CHOLESKY-S SQUARE-ROOT METHOD,
   WHICH GENERATES THE UPPER BAND MATRIX TU SUCH THAT
   TRANSPOSE(TU)*TU=A. TU IS RETURNED AS RESULT ON THE
   LOCATIONS OF A. FURTHER, DEPENDENT ON THE ACTUAL VALUE OF
   IOP, DIVISION OF R BY TRANSPOSE(TU) AND/OR TU IS PERFORMED
   AND THE RESULT IS RETURNED ON THE LOCATIONS OF R.
   FOR REFERENCE, SEE H. RUTISHAUSER, ALGORITHMUS 1 - LINEARES
   GLEICHUNGSSYSTEM MIT SYMMETRISCHER POSITIV-DEFINITER
   BANDMATRIX NACH CHOLESKY - , COMPUTING (ARCHIVES FOR
   ELECTRONIC COMPUTING), VOL.1, ISS.1 (1966), PP.77-78.