*
* If you have a subroutine for the matrix-vector product,
* this routine examplifies how to convert the matrix
* into the standard Harewell-Boeing format (compress column
* storage).
*
      INTEGER            NMAX, NZAMAX
      PARAMETER          ( NMAX = 500, NZAMAX = 60000 )
*
      DOUBLE PRECISION    ZERO, ONE
      PARAMETER           (ZERO = 0.0D+0, ONE = 1.0D+0 )
*
      CHARACTER          TITLE*72, KEY*8, TYPE*3
      INTEGER            IA(NZAMAX), JA(NMAX+1)
      DOUBLE PRECISION   A(NZAMAX), RHS(1)
      INTEGER            JOB, NROW, NCOL, K, IFMT, IOUNIT
*
      DOUBLE PRECISION    X(NMAX), Y(NMAX)
      INTEGER             N, I, J 
*
*     open a file to write the output matrix 
*
      open(unit=14,file='grcarmat',status='new')
*
      N = 100 
*
*     Harwell-Boeing format part
*
      IF( N .GT. NMAX )THEN
         PRINT*,'INCREASE NMAX TO ',N
         STOP
      ENDIF
*
      NROW = N
      NCOL = N
      K = 0
      JA(1) = 1
*
      DO 20 J = 1, N
*
*        compute column j of A
*
         DO 10 I = 1, N
            X(I) = ZERO
10       CONTINUE
         X(J) = ONE
*
         CALL MVMGRC(1,N,1,X,NMAX,Y,NMAX)
*
         DO 30 I = 1, N
            IF( Y(I) .NE. ZERO )THEN
               K = K + 1 
               A(K) = Y(i)
               IA(K) = I
            ENDIF
30       CONTINUE
         JA(J+1) = K + 1
*
20    CONTINUE
* 
      PRINT*,'DONE ALL MATVEC    N=  ',N,'  NNZ= ', JA(N+1)-1
*
      IF( JA(N+1) .GE. NZAMAX )THEN
         PRINT*,'INCREASE NZAMAX TO ',JA(N+1)-1
         STOP
      ENDIf
*
*     name the tile, key, type, etc 
*
      TITLE      ='GRCAR MATRIX'
      KEY        ='GRCAR100'
      TYPE       ='RUA'
      IFMT       = 16
      IOUNIT     = 14
      JOB        = 2
*
      CALL DPRTMT(NROW,NCOL,A,IA,JA,RHS,TITLE,KEY,TYPE,IFMT,JOB,IOUNIT)
*
      CLOSE( 14 )
*
      STOP
      end