#include "cppdefs.h" #ifdef NBQ !------------------------------------------------------------------------------ ! NHOMS ! Non Hydrostatic Ocean Modeling System !------------------------------------------------------------------------------ ! !> @note Main web documentation ! ! DESCRIPTION: ! !> @brief !! Sparskit library !! : Compressed Sparse Row format (CSR). !> @details Additional documentation: !! - Youcef Saad !! - Sparskit ! REVISION HISTORY: ! !> @authors !> @date 2015 January !> @todo ! !------------------------------------------------------------------------------ subroutine amubdg ( nrow, ncol, ncolb, ja, ia, jb, ib, ndegr, nnz, iw ) !*****************************************************************************80 ! !> AMUBDG gets the number of nonzero elements in each row of A * B. ! !> Discussion: !! !! The routine also computes the total number of nonzero elements in A * B. !! !! Method: A' * A = sum [over i = 1, nrow] a(i)^T a(i) !! where a(i) = i-th row of A. We must be careful not to add the !! elements already accounted for. ! !> Modified: !! !! 07 January 2004 ! !> Author: !! !! Youcef Saad ! !> Parameters: !! !! Input, integer NROW, the row dimension of the matrix A. !! !! Input, integer NCOL, the column dimension of the matrix A, !! (and the row dimension of B). !! !! Input, integer NCOLB, the column dimension of the matrix B. !! !! Input, ja, ia= row structure of input matrix A: ja = column indices of !! the nonzero elements of A stored by rows. !! ia = pointer to beginning of each row in ja. !! !! Input, jb, ib, the row structure of input matrix B: jb = column indices of !! the nonzero elements of A stored by rows. !! ib is a pointer to beginning of each row in jb. !! !! Output, integer NDEGR(NROW), contains the degrees (the number of !! nonzeros in each row of the matrix A * B. !! !! Output, integer NNZ, the number of nonzero elements found in A * B. !! !! Workspace, integer IW(NCOLB). !! implicit none integer ncol integer ncolb integer nrow integer ia(*) integer ib(*) integer ii integer iw(*) integer j integer ja(*) integer jb(*) integer jc integer jr integer k integer last integer ldg integer ndegr(*) integer nnz iw(1:ncolb) = 0 ndegr(1:nrow) = 0 do ii = 1, nrow ! ! For each row of A. ! ldg = 0 ! ! End-of-linked list. ! last = -1 do j = ia(ii), ia(ii+1)-1 ! ! Row number to be added. ! jr = ja(j) do k = ib(jr), ib(jr+1)-1 jc = jb(k) ! ! Add one element to the linked list. ! if ( iw(jc) == 0 ) then ldg = ldg + 1 iw(jc) = last last = jc end if end do end do ndegr(ii) = ldg ! ! Reset IW to zero. ! do k = 1, ldg j = iw(last) iw(last) = 0 last = j end do end do nnz = sum ( ndegr(1:nrow) ) return end subroutine amux ( n, x, y, a, ja, ia ) !subroutine amux ( n, x, y, a, ja, nl ) !*****************************************************************************80 ! !! AMUX multiplies a CSR matrix A times a vector. ! ! Discussion: ! ! This routine multiplies a matrix by a vector using the dot product form. ! Matrix A is stored in compressed sparse row storage. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), and array of length equal to the column dimension ! of A. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real Y(N), the product A * X. ! ! use module_parallele ! #MPI# use module_nh implicit none ! include 'mkl_spblas.h' integer n real ( kind = 8 ) a(*) integer i,j integer ia(*) integer ja(*) integer k ! integer nl ! character(len=1) :: TRANS='N' real ( kind = 8 ) t real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) do i = 1, n ! ! Compute the inner product of row I with vector X. ! t = 0.0D+00 do k = ia(i), ia(i+1)-1 ! do k = 1 + nl*(i-1), nl*i t = t + a(k) * x(ja(k)) end do y(i) = t end do return end subroutine amuxa ( n, x, y, a, ja, ia ) !subroutine amux ( n, x, y, a, ja, nl ) !*****************************************************************************80 ! !! AMUX multiplies a CSR matrix A times a vector. ! ! Discussion: ! ! This routine multiplies a matrix by a vector using the dot product form. ! Matrix A is stored in compressed sparse row storage. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), and array of length equal to the column dimension ! of A. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real Y(N), the product A * X. ! ! use module_parallele ! #MPI# use module_nh implicit none ! include 'mkl_spblas.h' integer n real ( kind = 8 ) a(*) integer i,j integer ia(*) integer ja(*) integer k ! integer nl ! character(len=1) :: TRANS='N' real ( kind = 8 ) t real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) do i = 1, n ! ! Compute the inner product of row I with vector X. ! t = 0.0D+00 do k = ia(i), ia(i+1)-1 ! do k = 1 + nl*(i-1), nl*i t = t + a(k) * x(ja(k)) end do y(i) = y(i)+t end do return end !*****************************************************************************80 ! amux pour MOM !*****************************************************************************80 subroutine amux_s ( n , n1, x, y, a, ja, ia ) ! use module_parallele ! #MPI# use module_nh implicit none ! include 'mkl_spblas.h' integer n,n1 real ( kind = 8 ) a(*) integer i integer ia(*) integer ja(*) integer k character(len=1) :: TRANS='N' real ( kind = 8 ) t real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) double precision :: cput1, cput2 !cput1=mpi_wtime() do i = 1, n1 t = 0.0D+00 do k = ia(i), ia(i+1)-1 t = t + a(k) * x(ja(k)) end do y(i) = t end do do i = n1+1, n y(i) = a(ia(i)) * x(ja(ia(i))) + a(ia(i)+1) * x(ja(ia(i)+1)) end do return end subroutine amub2 ( nrow, ncol, job, a, ja, ia, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !***************************************************************************** ! !! AMUB performs the matrix product C = A * B. ! ! Discussion: ! ! The column dimension of B is not needed. ! ! Modified: ! ! 08 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer NCOL, the column dimension of the matrix. ! ! Input, integer JOB, job indicator. When JOB = 0, only the structure ! is computed, that is, the arrays JC and IC, but the real values ! are ignored. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, b, jb, ib, matrix B in compressed sparse row format. ! ! Input, integer NZMAX, the length of the arrays c and jc. ! The routine will stop if the result matrix C has a number ! of elements that exceeds exceeds NZMAX. ! ! on return: ! ! c, ! jc, ! ic = resulting matrix C in compressed sparse row sparse format. ! ! ierr = integer. serving as error message. ! ierr = 0 means normal return, ! ierr > 0 means that amub stopped while computing the ! i-th row of C with i = ierr, because the number ! of elements in C exceeds nzmax. ! ! work arrays: ! ! iw = integer work array of length equal to the number of ! columns in A. ! implicit none integer ncol integer nrow integer nzmax real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(nzmax) ! integer ia(nrow+1) ! integer ib(ncol+1) ! integer ic(ncol+1) integer ia(*) integer ib(*) integer ic(*) integer ierr integer ii integer iw(*) integer ja(*) integer jb(*) integer jc(nzmax) integer jcol integer jj integer job integer jpos integer k integer ka integer kb integer len real ( kind = 8 ) scal logical values values = ( job /= 0 ) len = 0 ic(1) = 1 ierr = 0 ! ! Initialize IW. ! iw(1:ncol) = 0 do ii = 1, nrow ! ! Row I. ! do ka = ia(ii), ia(ii+1)-1 if ( values ) then scal = a(ka) end if jj = ja(ka) do kb = ib(jj), ib(jj+1)-1 jcol = jb(kb) jpos = iw(jcol) if ( jpos == 0 ) then len = len + 1 if ( nzmax < len ) then ierr = ii return end if jc(len) = jcol iw(jcol)= len if ( values ) then c(len) = scal * b(kb) end if else if ( values ) then c(jpos) = c(jpos) + scal * b(kb) end if end if end do end do do k = ic(ii), len iw(jc(k)) = 0 end do ic(ii+1) = len + 1 end do return end subroutine amub3 ( nrow, ncol, job, alpha, a, ja, ia, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !***************************************************************************** ! !! AMUB performs the matrix product C = A * B. ! ! Discussion: ! ! The column dimension of B is not needed. ! ! Modified: ! ! 08 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer NCOL, the column dimension of the matrix. ! ! Input, integer JOB, job indicator. When JOB = 0, only the structure ! is computed, that is, the arrays JC and IC, but the real values ! are ignored. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, b, jb, ib, matrix B in compressed sparse row format. ! ! Input, integer NZMAX, the length of the arrays c and jc. ! The routine will stop if the result matrix C has a number ! of elements that exceeds exceeds NZMAX. ! ! on return: ! ! c, ! jc, ! ic = resulting matrix C in compressed sparse row sparse format. ! ! ierr = integer. serving as error message. ! ierr = 0 means normal return, ! ierr > 0 means that amub stopped while computing the ! i-th row of C with i = ierr, because the number ! of elements in C exceeds nzmax. ! ! work arrays: ! ! iw = integer work array of length equal to the number of ! columns in A. ! implicit none integer ncol integer nrow integer nzmax real ( kind = 8 ) alpha real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(nzmax) ! integer ia(nrow+1) ! integer ib(ncol+1) ! integer ic(ncol+1) integer ia(*) integer ib(*) integer ic(*) integer ierr integer ii integer iw(*) integer ja(*) integer jb(*) integer jc(nzmax) integer jcol integer jj integer job integer jpos integer k integer ka integer kb integer len real ( kind = 8 ) scal logical values values = ( job /= 0 ) len = 0 ic(1) = 1 ierr = 0 ! ! Initialize IW. ! iw(1:ncol) = 0 do ii = 1, nrow ! ! Row I. ! do ka = ia(ii), ia(ii+1)-1 if ( values ) then scal = a(ka)*alpha end if jj = ja(ka) do kb = ib(jj), ib(jj+1)-1 jcol = jb(kb) jpos = iw(jcol) if ( jpos == 0 ) then len = len + 1 if ( nzmax < len ) then ierr = ii return end if jc(len) = jcol iw(jcol)= len if ( values ) then c(len) = scal * b(kb) end if else c(jpos) = c(jpos) + scal * b(kb) end if end do end do do k = ic(ii), len iw(jc(k)) = 0 end do ic(ii+1) = len + 1 end do return end subroutine aplb ( nrow, ncol, job, a, ja, ia, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !*****************************************************************************80 ! !! APLB performs the CSR matrix sum C = A + B. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer ( kind = 4 ) NROW, the row dimension of A and B. ! ! Input, integer ( kind = 4 ) NCOL, the column dimension of A and B. ! ! Input, integer ( kind = 4 ) JOB. When JOB = 0, only the structure ! (i.e. the arrays jc, ic) is computed and the ! real values are ignored. ! ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format. ! ! nzmax = integer ( kind = 4 ). The length of the arrays c and jc. ! amub will stop if the result matrix C has a number ! of elements that exceeds exceeds nzmax. See ierr. ! ! on return: ! ! c, ! jc, ! ic = resulting matrix C in compressed sparse row sparse format. ! ! ierr = integer ( kind = 4 ). serving as error message. ! ierr = 0 means normal return, ! ierr > 0 means that amub stopped while computing the ! i-th row of C with i = ierr, because the number ! of elements in C exceeds nzmax. ! ! work arrays: ! ! iw = integer ( kind = 4 ) work array of length equal to the number of ! columns in A. ! implicit none integer ( kind = 4 ) ncol integer ( kind = 4 ) nrow real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer ( kind = 4 ) ia(nrow+1) integer ( kind = 4 ) ib(nrow+1) integer ( kind = 4 ) ic(nrow+1) integer ( kind = 4 ) ierr integer ( kind = 4 ) ii integer ( kind = 4 ) iw(ncol) integer ( kind = 4 ) ja(*) integer ( kind = 4 ) jb(*) integer ( kind = 4 ) jc(*) integer ( kind = 4 ) jcol integer ( kind = 4 ) job integer ( kind = 4 ) jpos integer ( kind = 4 ) k integer ( kind = 4 ) ka integer ( kind = 4 ) kb integer ( kind = 4 ) len integer ( kind = 4 ) nzmax logical values values = ( job /= 0 ) ierr = 0 len = 0 ic(1) = 1 iw(1:ncol) = 0 do ii = 1, nrow ! ! Row I. ! do ka = ia(ii), ia(ii+1)-1 len = len + 1 jcol = ja(ka) if ( nzmax < len ) then ierr = ii return end if jc(len) = jcol if ( values ) then c(len) = a(ka) end if iw(jcol) = len end do do kb = ib(ii), ib(ii+1)-1 jcol = jb(kb) jpos = iw(jcol) if ( jpos == 0 ) then len = len + 1 if ( nzmax < len ) then ierr = ii return end if jc(len) = jcol if ( values ) then c(len) = b(kb) end if iw(jcol)= len else if ( values ) then c(jpos) = c(jpos) + b(kb) end if end if end do do k = ic(ii), len iw(jc(k)) = 0 end do ic(ii+1) = len+1 end do return end subroutine aplb1 ( nrow, ncol, job, a, ja, ia, b, jb, ib, c, jc, ic, & nzmax, ierr ) !*****************************************************************************80 ! !! APLB1 performs the sum C = A + B for sorted CSR matrices. ! ! Discussion: ! ! The difference between this routine and APLB is that here the ! resulting matrix is such that the elements of each row are sorted, ! with increasing column indices in each row, provided the original ! matrices are sorted in the same way. ! ! This routine will not work if either of the two input matrices is ! not sorted. ! ! Modified: ! ! 11 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer ( kind = 4 ) NROW, the row dimension of A and B. ! ! Input, integer ( kind = 4 ) NCOL, the column dimension of A and B. ! ! Input, integer ( kind = 4 ) JOB. When JOB = 0, only the structure ! (i.e. the arrays jc, ic) is computed and the ! real values are ignored. ! ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format with entries sorted. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format with entries sorted ! ascendly in each row ! ! nzmax = integer ( kind = 4 ). The length of the arrays c and jc. ! amub will stop if the result matrix C has a number ! of elements that exceeds exceeds nzmax. See ierr. ! ! on return: ! ! c, ! jc, ! ic = resulting matrix C in compressed sparse row sparse format ! with entries sorted ascendly in each row. ! ! ierr = integer ( kind = 4 ). serving as error message. ! ierr = 0 means normal return, ! ierr > 0 means that amub stopped while computing the ! i-th row of C with i = ierr, because the number ! of elements in C exceeds nzmax. ! implicit none integer ( kind = 4 ) nrow real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer ( kind = 4 ) i integer ( kind = 4 ) ia(nrow+1) integer ( kind = 4 ) ib(nrow+1) integer ( kind = 4 ) ic(nrow+1) integer ( kind = 4 ) ierr integer ( kind = 4 ) j1 integer ( kind = 4 ) j2 integer ( kind = 4 ) ja(*) integer ( kind = 4 ) jb(*) integer ( kind = 4 ) jc(*) integer ( kind = 4 ) job integer ( kind = 4 ) ka integer ( kind = 4 ) kamax integer ( kind = 4 ) kb integer ( kind = 4 ) kbmax integer ( kind = 4 ) kc integer ( kind = 4 ) ncol integer ( kind = 4 ) nzmax logical values values = ( job /= 0 ) ierr = 0 kc = 1 ic(1) = kc do i = 1, nrow ka = ia(i) kb = ib(i) kamax = ia(i+1) - 1 kbmax = ib(i+1) - 1 do if ( ka <= kamax ) then j1 = ja(ka) else j1 = ncol + 1 end if if ( kb <= kbmax ) then j2 = jb(kb) else j2 = ncol + 1 end if ! ! Three cases ! if ( j1 == j2 ) then if ( values ) then c(kc) = a(ka) + b(kb) end if jc(kc) = j1 ka = ka + 1 kb = kb + 1 kc = kc + 1 else if ( j1 < j2 ) then jc(kc) = j1 if ( values ) then c(kc) = a(ka) end if ka = ka + 1 kc = kc + 1 else if ( j2 < j1 ) then jc(kc) = j2 if ( values ) then c(kc) = b(kb) end if kb = kb + 1 kc = kc + 1 end if if ( nzmax < kc ) then ierr = i return end if if ( kamax < ka .and. kbmax < kb ) then exit end if end do ic(i+1) = kc end do return end subroutine aplsca1 ( nrow, a, ja, ia, iw ) !*****************************************************************************80 ! !! APLSCA adds 1. to the diagonal entries of a sparse matrix A :=A + s I ! ! Discussion: ! ! The column dimension of A is not needed. ! ! important: the matrix A may be expanded slightly to allow for ! additions of nonzero elements to previously nonexisting diagonals. ! The is no checking as to whether there is enough space appended ! to the arrays a and ja. if not sure allow for n additional ! elements. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer ( kind = 4 ) NROW, the row dimension of the matrix. ! ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! on return: ! ! ! a, ! ja, ! ia = matrix A with diagonal elements shifted (or created). ! ! iw = integer ( kind = 4 ) work array of length n. On return iw will ! contain the positions of the diagonal entries in the ! output matrix. (i.e., a(iw(k)), ja(iw(k)), k = 1,...n, ! are the values/column indices of the diagonal elements ! of the output matrix. ). ! implicit none integer ( kind = 4 ) nrow real ( kind = 8 ) a(*) integer ( kind = 4 ) ia(nrow+1) integer ( kind = 4 ) icount integer ( kind = 4 ) ii integer ( kind = 4 ) iw(*) integer ( kind = 4 ) j integer ( kind = 4 ) ja(*) integer ( kind = 4 ) k integer ( kind = 4 ) k1 integer ( kind = 4 ) k2 integer ( kind = 4 ) ko logical test call diapos ( nrow, ja, ia, iw ) icount = 0 do j = 1, nrow if ( iw(j) == 0 ) then icount = icount + 1 else a(iw(j)) = a(iw(j)) + 1. end if end do ! ! If no diagonal elements to insert in data structure, return. ! if ( icount == 0 ) then return end if ! ! Shift the nonzero elements if needed, to allow for created ! diagonal elements. ! ko = ia(nrow+1) + icount ! ! Copy rows backward. ! do ii = nrow, 1, -1 ! ! Go through row II. ! k1 = ia(ii) k2 = ia(ii+1) - 1 ia(ii+1) = ko test = ( iw(ii) == 0 ) do k = k2, k1, -1 j = ja(k) if ( test .and. j < ii ) then test = .false. ko = ko - 1 a(ko) = 1. ja(ko) = ii iw(ii) = ko end if ko = ko - 1 a(ko) = a(k) ja(ko) = j end do ! ! The diagonal element has not been added yet. ! if ( test ) then ko = ko - 1 a(ko) = 1. ja(ko) = ii iw(ii) = ko end if end do ia(1) = ko return end subroutine diapos ( n, ja, ia, idiag ) !*****************************************************************************80 ! !! DIAPOS returns the positions of the diagonal elements of a sparse matrix. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the row dimension of the matrix. ! ! Input, JA(*), IA(N+1), the matrix information, (but no values) ! in CSR Compressed Sparse Row format. ! ! Output, integer ( kind = 4 ) IDIAG(N); the I-th entry of IDIAG points to the ! diagonal element A(I,I) in the arrays A and JA. That is, ! A(IDIAG(I)) = element A(I,I) of matrix A. If no diagonal element ! is found, the entry is set to 0. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) ia(n+1) integer ( kind = 4 ) idiag(n) integer ( kind = 4 ) ja(*) integer ( kind = 4 ) k idiag(1:n) = 0 ! ! Sweep through the data structure. ! do i = 1, n do k = ia(i), ia(i+1) -1 if ( ja(k) == i ) then idiag(i) = k end if end do end do return end #ifdef FOR_TESTS subroutine amux_1 ( n, x, y, a, ja, ia ) !subroutine amux ( n, x, y, a, ja, nl ) !*****************************************************************************80 ! !! AMUX multiplies a CSR matrix A times a vector. ! ! Discussion: ! ! This routine multiplies a matrix by a vector using the dot product form. ! Matrix A is stored in compressed sparse row storage. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), and array of length equal to the column dimension ! of A. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real Y(N), the product A * X. ! ! use module_parallele ! #MPI# use module_nh use module_nbq implicit none ! include 'mkl_spblas.h' integer n real ( kind = 8 ) a(*) integer i,j integer ia(*) integer ja(*) integer k ! integer nl ! character(len=1) :: TRANS='N' real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) namux1_nbq=0 do i = 1, n ! ! Compute the inner product of row I with vector X. ! ! t = 0.0D+00 do k = ia(i), ia(i+1)-1 ! do k = 1 + nl*(i-1), nl*i ! t = t + a(k) * x(ja(k)) namux1_nbq = namux1_nbq+1 amux1a_nbq(namux1_nbq)=i amux1b_nbq(namux1_nbq)=k amux1c_nbq(namux1_nbq)=ja(k) end do ! y(i) = t end do return end subroutine amux_2 ( n, x, y, a ) !subroutine amux ( n, x, y, a, ja, nl ) !*****************************************************************************80 ! !! AMUX multiplies a CSR matrix A times a vector. ! ! Discussion: ! ! This routine multiplies a matrix by a vector using the dot product form. ! Matrix A is stored in compressed sparse row storage. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), and array of length equal to the column dimension ! of A. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real Y(N), the product A * X. ! ! use module_parallele ! #MPI# use module_nh use module_nbq implicit none ! include 'mkl_spblas.h' integer n real ( kind = 8 ) a(*) integer i,j integer k ! integer nl ! character(len=1) :: TRANS='N' real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) integer icall y(1:n)=0. do l_nbq=1,namux1_nbq y(amux1a_nbq(l_nbq)) = y(amux1a_nbq(l_nbq)) +a(amux1b_nbq(l_nbq))*x(amux1c_nbq(l_nbq)) enddo return end #endif #endif