*DECK CPOIR SUBROUTINE CPOIR (A, LDA, N, V, ITASK, IND, WORK) C***BEGIN PROLOGUE CPOIR C***PURPOSE Solve a positive definite Hermitian system of linear C equations. Iterative refinement is used to obtain an C error estimate. C***LIBRARY SLATEC C***CATEGORY D2D1B C***TYPE COMPLEX (SPOIR-S, CPOIR-C) C***KEYWORDS HERMITIAN, LINEAR EQUATIONS, POSITIVE DEFINITE, SYMMETRIC C***AUTHOR Voorhees, E. A., (LANL) C***DESCRIPTION C C Subroutine CPOIR solves a complex positive definite Hermitian C NxN system of single precision linear equations using LINPACK C subroutines CPOFA and CPOSL. One pass of iterative refine- C ment is used only to obtain an estimate of the accuracy. That C is, if A is an NxN complex positive definite Hermitian matrix C and if X and B are complex N-vectors, then CPOIR solves the C equation C C A*X=B. C C Care should be taken not to use CPOIR with a non-Hermitian C matrix. C C The matrix A is first factored into upper and lower C triangular matrices R and R-TRANSPOSE. These C factors are used to calculate the solution, X. C Then the residual vector is found and used C to calculate an estimate of the relative error, IND. C IND estimates the accuracy of the solution only when the C input matrix and the right hand side are represented C exactly in the computer and does not take into account C any errors in the input data. C C If the equation A*X=B is to be solved for more than one vector C B, the factoring of A does not need to be performed again and C the option to only solve (ITASK .GT. 1) will be faster for C the succeeding solutions. In this case, the contents of A, C LDA, N, and WORK must not have been altered by the user C following factorization (ITASK=1). IND will not be changed C by CPOIR in this case. C C Argument Description *** C A COMPLEX(LDA,N) C the doubly subscripted array with dimension (LDA,N) C which contains the coefficient matrix. Only the C upper triangle, including the diagonal, of the C coefficient matrix need be entered. A is not C altered by the routine. C LDA INTEGER C the leading dimension of the array A. LDA must be great- C er than or equal to N. (terminal error message IND=-1) C N INTEGER C the order of the matrix A. N must be greater than C or equal to one. (terminal error message IND=-2) C V COMPLEX(N) C on entry, the singly subscripted array(vector) of di- C mension N which contains the right hand side B of a C system of simultaneous linear equations A*X=B. C on return, V contains the solution vector, X . C ITASK INTEGER C if ITASK = 1, the matrix A is factored and then the C linear equation is solved. C if ITASK .GT. 1, the equation is solved using the existing C factored matrix A (stored in WORK). C if ITASK .LT. 1, then terminal terminal error IND=-3 is C printed. C IND INTEGER C GT. 0 IND is a rough estimate of the number of digits C of accuracy in the solution, X. IND=75 means C that the solution vector X is zero. C LT. 0 see error message corresponding to IND below. C WORK COMPLEX(N*(N+1)) C a singly subscripted array of dimension at least N*(N+1). C C Error Messages Printed *** C C IND=-1 terminal N is greater than LDA. C IND=-2 terminal N is less than one. C IND=-3 terminal ITASK is less than one. C IND=-4 terminal The matrix A is computationally singular C or is not positive definite. C A solution has not been computed. C IND=-10 warning The solution has no apparent significance. C the solution may be inaccurate or the matrix C a may be poorly scaled. C C NOTE- the above terminal(*fatal*) error messages are C designed to be handled by XERMSG in which C LEVEL=1 (recoverable) and IFLAG=2 . LEVEL=0 C for warning error messages from XERMSG. Unless C the user provides otherwise, an error message C will be printed followed by an abort. C C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. C Stewart, LINPACK Users' Guide, SIAM, 1979. C***ROUTINES CALLED CCOPY, CPOFA, CPOSL, DCDOT, R1MACH, SCASUM, XERMSG C***REVISION HISTORY (YYMMDD) C 800530 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900510 Convert XERRWV calls to XERMSG calls, cvt GOTO's to C IF-THEN-ELSE. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE CPOIR C INTEGER LDA,N,ITASK,IND,INFO,J COMPLEX A(LDA,*),V(*),WORK(N,*) REAL SCASUM,XNORM,DNORM,R1MACH DOUBLE PRECISION DR1,DI1,DR2,DI2 CHARACTER*8 XERN1, XERN2 C***FIRST EXECUTABLE STATEMENT CPOIR IF (LDA.LT.N) THEN IND = -1 WRITE (XERN1, '(I8)') LDA WRITE (XERN2, '(I8)') N CALL XERMSG ('SLATEC', 'CPOIR', 'LDA = ' // XERN1 // * ' IS LESS THAN N = ' // XERN2, -1, 1) RETURN ENDIF C IF (N.LE.0) THEN IND = -2 WRITE (XERN1, '(I8)') N CALL XERMSG ('SLATEC', 'CPOIR', 'N = ' // XERN1 // * ' IS LESS THAN 1', -2, 1) RETURN ENDIF C IF (ITASK.LT.1) THEN IND = -3 WRITE (XERN1, '(I8)') ITASK CALL XERMSG ('SLATEC', 'CPOIR', 'ITASK = ' // XERN1 // * ' IS LESS THAN 1', -3, 1) RETURN ENDIF C IF (ITASK.EQ.1) THEN C C MOVE MATRIX A TO WORK C DO 10 J=1,N CALL CCOPY(N,A(1,J),1,WORK(1,J),1) 10 CONTINUE C C FACTOR MATRIX A INTO R C CALL CPOFA(WORK,N,N,INFO) C C CHECK FOR SINGULAR OR NOT POS.DEF. MATRIX C IF (INFO.NE.0) THEN IND = -4 CALL XERMSG ('SLATEC', 'CPOIR', * 'SINGULAR OR NOT POSITIVE DEFINITE - NO SOLUTION', -4, 1) RETURN ENDIF ENDIF C C SOLVE AFTER FACTORING C MOVE VECTOR B TO WORK C CALL CCOPY(N,V(1),1,WORK(1,N+1),1) CALL CPOSL(WORK,N,N,V) C C FORM NORM OF X0 C XNORM = SCASUM(N,V(1),1) IF (XNORM.EQ.0.0) THEN IND = 75 RETURN ENDIF C C COMPUTE RESIDUAL C DO 40 J=1,N CALL DCDOT(J-1,-1.D0,A(1,J),1,V(1),1,DR1,DI1) CALL DCDOT(N-J+1,1.D0,A(J,J),LDA,V(J),1,DR2,DI2) DR1 = DR1+DR2-DBLE(REAL(WORK(J,N+1))) DI1 = DI1+DI2-DBLE(AIMAG(WORK(J,N+1))) WORK(J,N+1) = CMPLX(REAL(DR1),REAL(DI1)) 40 CONTINUE C C SOLVE A*DELTA=R C CALL CPOSL(WORK,N,N,WORK(1,N+1)) C C FORM NORM OF DELTA C DNORM = SCASUM(N,WORK(1,N+1),1) C C COMPUTE IND (ESTIMATE OF NO. OF SIGNIFICANT DIGITS) C AND CHECK FOR IND GREATER THAN ZERO C IND = -LOG10(MAX(R1MACH(4),DNORM/XNORM)) IF (IND.LE.0) THEN IND = -10 CALL XERMSG ('SLATEC', 'CPOIR', * 'SOLUTION MAY HAVE NO SIGNIFICANCE', -10, 0) ENDIF RETURN END