getmatvec.f

网络带宽测试工具

      Subroutine getmatvec( n, nelmts, indx, rowp, matvals, invec )
! ----------------------------------------------------------------------
! --- 'getmatvec' makes the data of the CRS matrix and the input vector.
!     The number of non-zero elements/row may randomly vary by 20%.
! ----------------------------------------------------------------------
      Use                      numerics
      Use                      ran_module
      Implicit                 None

      Integer               :: n, nelmts
      Integer               :: indx(nelmts), rowp(n+1)
      Real(l_)              :: matvals(nelmts), invec(n)

      Real(l_)              :: varelts(nelmts)
      Integer               :: i, j, nelm
      Integer               :: offset, shindx
      Real(l_)              :: epr, offdiag
! ----------------------------------------------------------------------
! --- Generate the rowpointers.
! ----------------------------------------------------------------------
      x1 = 2006; x2 = 2003                        !<-- Set random seeds.
      Call rinit
      epr = Real( nelmts, l_ )/Real(n, l_ )
      rowp(1) = 1
      Do i = 2, n
         rowp(i) = rowp(i-1) +
     &             Ceiling( epr*( 1.0 - 0.2*dran0() ) )
      End Do
      rowp(n+1) = nelmts
! ----------------------------------------------------------------------
! --- Generate the index array (and make sure there is a diag. element).
! ----------------------------------------------------------------------     
      Call ranfil( varelts, nelmts )
      indx = n*varelts + 1
      Do i = 1, n + 1
         indx(rowp(i)) = i
      End Do
! ----------------------------------------------------------------------
! --- Sort indices per row and adjust multiple indentical row entries.
! ----------------------------------------------------------------------     
      Do i = 1, n
         Call iqsort( indx, nelmts, rowp(i), rowp(i+1) - 1 )
         nelm = rowp(i+1) - rowp(i)
         offdiag = 1.0_l_/Real( Max( nelm - 1, 1 ), l_ )
         matvals(rowp(i):rowp(i+1)-1) = -0.97_l_*offdiag
         Do j = rowp(i), rowp(i+1) - 1
            If ( indx(j) == indx(j+1) ) Then
               If ( indx(j+1) >= n ) Then
                  indx(j+1) = indx(j+1) - 2
               Else
                  indx(j+1) = indx(j+1) + 1
               End If
            End If
         End Do
! ----------------------------------------------------------------------
! --- Set diagonal elements to 1. The RHS is the sum of row elements.
!     So, x(i) should converge to 1 for row(i), (i = 1, n).
! ----------------------------------------------------------------------
         Call diagsk( nelm, indx(rowp(i)), i, offset )
         shindx = rowp(i) + offset
         indx(shindx) = i
         matvals(shindx) = 1.0_l_
         invec(i) = Sum( matvals(rowp(i):rowp(i+1)-1) )
      End Do
! ----------------------------------------------------------------------
      End Subroutine getmatvec