      program example
*
*     example program solving ax=b via scalapack routine pdgesv
*

      integer n, nb, nprow, npcol      
      parameter (n = 4000, nb = 64, nprow =2, npcol = 3)

      integer nlocr, nlocc
      parameter (nlocr = 2048, nlocc=1344)

      integer            mycol, myrow
      integer            ictxt, info
      double precision   anorm, bnorm, eps, resid, xnorm

      integer            desca(9), descb(9)
      double precision   a(nlocr, nlocc), a0(nlocr, nlocc)
      double precision   b(nlocr, 1), b0(nlocr, 1)
      double precision   work(nlocr)
      integer            ipiv(nlocr+nb)


*****

*     initialize the process grid
      call sl_init( ictxt, nprow, npcol )
      call blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )

*     initialize the array descriptors for the matrices a and b
      call descinit( desca, n, n, nb, nb, 0, 0, ictxt, nlocr, info )
      call descinit( descb, n, 1, nb, 1, 0, 0, ictxt, nlocr, info )

*     generate distributed matrices a and b
      call matinit( a, desca, b, descb )

*     make a copy of a and b for checking purposes
      call pdlacpy( 'a', n, n, a, 1, 1, desca, a0, 1, 1, desca )
      call pdlacpy( 'a', n, 1, b, 1, 1, descb, b0, 1, 1, descb )

*     solve the linear system a * x = b
      call pdgesv( n, 1, a, 1, 1, desca, ipiv, b, 1, 1, descb, info )

*****

*     compute residual ||a*x  - b||/( ||x|| * ||a|| * eps * n )

*     get machine precision
      eps = pdlamch( ictxt, 'epsilon' )

*     get infinity-norm of a and b
      anorm = pdlange( 'i', n, n, a, 1, 1, desca, work )
      bnorm = pdlange( 'i', n, 1, b, 1, 1, descb, work )

*     compute a0*b - b0  (i.e. a*x - b) and its norm
      call pdgemm( 'n', 'n', n, 1, n, 1.0d0, a0, 1, 1, desca, b, 1, 1,
     $             descb, -1.0d0, b0, 1, 1, descb )
      xnorm = pdlange( 'i', n, 1, b0, 1, 1, descb, work )

*     check residual error
      resid = xnorm / ( anorm*bnorm*eps*dble( n ) )
      if( myrow.eq.0 .and. mycol.eq.0 ) then
         if( resid .lt. 10.0d+0 ) then
            print *, 'solution is correct! residual:', resid
         else
            print *, 'solution is incorrect! residual:', resid
         end if
      end if

      call blacs_exit( 0 )

      stop
      end

