*
* $Id: f010r.F,v 1.1.1.1 1996/02/15 17:48:42 mclareni Exp $
*
* $Log: f010r.F,v $
* Revision 1.1.1.1  1996/02/15 17:48:42  mclareni
* Kernlib
*
*
#include "kernnumt/pilot.h"
      SUBROUTINE F010R(N,K,KMAX,IDIM,A,AREF,B,BREF,W,R,RESID)
      REAL R(N)
      REAL U,S,A,AREF,B,BREF,W
      DIMENSION A(IDIM,N),AREF(IDIM,N),B(IDIM,2),BREF(IDIM,2),W(N)
C
C  REAL VERSION OF TEST FOR F010.
C  SETS RESID EQUAL TO THE LARGEST RESIDUAL FOUND WHEN INVERTING A
C  MATRIX OR SOLVING A LINEAR SYSTEM.
C
C     N        ORDER OF COEFFICIENT MATRIX.
C     K        CURRENT NUMBER OF RIGHT-HAND SIDES.
C     KMAX     MAXIMUM NUMBER OF RIGHT-HAND SIDES.
C     IDIM     FIRST DIMENSION OF THE 2D ARRAYS.
C     A,AREF   (REAL) 2D ARRAYS WITH AT LEAST N COLUMNS.
C     B,BREF   (REAL) 2D ARRAYS WITH AT LEAST KMAX COLUMNS.
C     W        (REAL) 1D ARRAY WITH AT LEAST N COLUMNS.
C     R        (REAL) 1D ARRAY WITH AT LEAST N COLUMNS.
C     RESID    (REAL) OUTPUT PARAMETER.
C
C  CALLS ... ALL REAL ENTRIES IN F010.
C        ... CERN PACKAGES F002 AND F003.
C        ... TEST-ROUTINE F010MR.
C
C  START.
      IF(K.EQ.0) GO TO 4
C
C  SET COPIES OF A RANDOM DIAGONALLY-DOMINANT SQUARE MATRIX IN A AND
C  AREF.
      U=1.
      S=FLOAT(N)
      CALL RMRAN(N,N,-U,U,A,A(1,2),A(2,1))
      CALL RVSET(N,S,W(1),W(2))
      CALL RVADD(N,A(1,1),A(2,2),W(1),W(2),A(1,1),A(2,2))
      CALL RMCPY(N,N,A,A(1,2),A(2,1),AREF,AREF(1,2),AREF(2,1))
C
C  SET COPIES OF A RANDOM RIGHT-HAND SIDE MATRIX IN B AND BREF.
      CALL RMRAN(N,KMAX,-U,U,B,B(1,2),B(2,1))
      CALL RMCPY(N,KMAX,B,B(1,2),B(2,1),BREF,BREF(1,2),BREF(2,1))
C
C  CALL F010 TO REPLACE A BY AINV AND B BY AINV*B. COMPUTE RESIDUALS.
      RESID=0.
      DO 3 ISUB=1,3
         CALL RMCPY(N,N,AREF,AREF(1,2),AREF(2,1),A,A(1,2),A(2,1))
         CALL RMCPY(N,KMAX,BREF,BREF(1,2),BREF(2,1),B,B(1,2),B(2,1))
         IF(ISUB.EQ.1) CALL REQINV(N,A,IDIM,R,IFAIL,K,B)
         IF(ISUB.EQ.2) CALL   RINV(N,A,IDIM,R,IFAIL)
         IF(ISUB.EQ.3) CALL   REQN(N,A,IDIM,R,IFAIL,K,B)
         RESINV=0.
         RESEQN=0.
         IF(ISUB.EQ.3) GO TO 1
C           (MAXIMUM RESIDUAL FOR INVERSION.)
            CALL RMMLT(N,N,N,AREF,AREF(1,2),AREF(2,1),
     *                   A,A(1,2),A(2,1),A,A(1,2),A(2,1),W)
            S=1.
            CALL RVSET(N,S,W(1),W(2))
            CALL RVSUB(N,A(1,1),A(2,2),W(1),W(2),A(1,1),A(2,2))
            RESINV=F010MR(N,N,IDIM,A)
    1    IF(ISUB.EQ.2) GO TO 2
C           (MAXIMUM RESIDUAL FOR LINEAR SYSTEM.)
            IF(K.NE.0) CALL RMMLT(N,N,K,AREF,AREF(1,2),AREF(2,1),
     *                              B,B(1,2),B(2,1),B,B(1,2),B(2,1),W)
            CALL RMSUB(N,KMAX,B,B(1,2),B(2,1),
     *                   BREF,BREF(1,2),BREF(2,1),B,B(1,2),B(2,1))
            RESEQN=F010MR(N,K,IDIM,B)
C        (FINAL RESIDUAL.)
    2    RESID=AMAX1(RESID,RESINV,RESEQN)
         IF(IFAIL.NE.0) RESID=99999.
    3 CONTINUE
      RETURN
C
C  TEST IFAIL FLAG.
    4 S=1.
      CALL RMSET(2,2,S,A,A(1,2),A(2,1))
      CALL RINV(2,A,IDIM,R,IFAIL)
      RESID=0.
      IF(IFAIL.NE.-1) RESID=99999.
      RETURN
      END