subroutine amask ( nrow, ncol, a, ja, ia, jmask, imask, c, jc, ic, iw, & nzmax, ierr ) !*****************************************************************************80 ! !! AMASK extracts a sparse matrix from a masked input matrix. ! ! Discussion: ! ! The routine looks at the positions defined by MASK, JMASK and IMASK. ! ! The algorithm is "in place": C, JC, IC can be the same as ! A, JA, IA. ! ! Modified: ! ! 08 January 2004 ! ! Author: ! ! Youcef Saad ! ! Reference: ! ! Youcef Saad, ! Sparsekit: a basic tool kit for sparse matrix computations, ! Technical Report, Computer Science Department, ! University of Minnesota, June 1994 ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer NCOL, the column dimension of the matrix. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, integer JMASK(*), IMASK((NROW+1), defining mask (pattern only) ! stored in compressed sparse row format. ! ! Input, integer NZMAX, the length of arrays C and JC. ! ! Output, c, jc, ic, the output matrix in Compressed Sparse Row format. ! ! Workspace, logical IW(NCOL). ! ! Input, integer NZMAX, the dimension of C. ! ! Output, integer IERR, serving as error message. ! ierr = 1 means normal return ! ierr > 1 means that amask stopped when processing ! row number ierr, because there was not enough space in ! c, jc according to the value of nzmax. ! implicit none integer ncol integer nrow integer nzmax real ( kind = 8 ) a(*) real ( kind = 8 ) c(nzmax) integer ia(nrow+1) integer ic(nrow+1) integer ierr integer ii integer imask(nrow+1) logical iw(ncol) integer j integer ja(*) integer jc(nzmax) integer jmask(*) integer k integer k1 integer k2 integer len ierr = 0 len = 0 iw(1:ncol) = .false. ! ! Unpack the mask for row II in IW. ! do ii = 1, nrow ! ! Save pointer in order to be able to do things in place. ! do k = imask(ii), imask(ii+1)-1 iw(jmask(k)) = .true. end do ! ! Add unmasked elemnts of row II. ! k1 = ia(ii) k2 = ia(ii+1)-1 ic(ii) = len+1 do k = k1, k2 j = ja(k) if ( iw(j) ) then len = len + 1 if ( nzmax < len ) then ierr = ii return end if jc(len) = j c(len) = a(k) end if end do do k = imask(ii), imask(ii+1)-1 iw(jmask(k)) = .false. end do end do ic(nrow+1) = len + 1 return end subroutine amub ( nrow, ncol, job, a, ja, ia, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !*****************************************************************************80 ! !! 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 ierr integer ii integer iw(ncol) 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 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(nrow+1) integer ib(ncol+1) integer ii integer iw(ncolb) integer j integer ja(*) integer jb(*) integer jc integer jr integer k integer last integer ldg integer ndegr(nrow) 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 amudia ( nrow, job, a, ja, ia, diag, b, jb, ib ) !*****************************************************************************80 ! !! AMUDIA performs the matrix by matrix product B = A * Diag (in place) ! ! Discussion: ! ! The column dimension of A is not needed. ! The algorithm is "in place", so B can take the place of A. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer JOB, job indicator. Job=0 means get array b only ! job = 1 means get b, and the integer arrays ib, jb. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, real DIAG(NROW), the diagonal matrix stored as a vector. ! ! Output, B(*), JB(*), IB(NROW+1), the resulting matrix B in ! compressed sparse row sparse format. ! implicit none integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) diag(nrow) integer ia(nrow+1) integer ib(nrow+1) integer ii integer ja(*) integer jb(*) integer job integer k integer k1 integer k2 do ii = 1, nrow ! ! Scale each element. ! k1 = ia(ii) k2 = ia(ii+1) - 1 do k = k1, k2 b(k) = a(k) * diag(ja(k)) end do end do if ( job == 0 ) then return end if ib(1) = ia(1) do ii = 1, nrow ib(ii) = ia(ii) do k = ia(ii), ia(ii+1)-1 jb(k) = ja(k) end do end do return end subroutine amux ( n, x, y, a, ja, ia ) !*****************************************************************************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. ! implicit none integer n real ( kind = 8 ) a(*) integer i integer ia(*) integer ja(*) integer k real ( kind = 8 ) t real ( kind = 8 ) x(*) real ( kind = 8 ) y(n) 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 t = t + a(k) * x(ja(k)) end do y(i) = t end do return end subroutine amuxd ( n, x, y, diag, ndiag, idiag, ioff ) !*****************************************************************************80 ! !! AMUXD multiplies a DIA matrix times a vector. ! ! Discussion: ! ! This routine multiplies a matrix by a vector when the original matrix ! is stored in the DIA diagonal storage format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), array of length equal to the column dimension of ! the A matrix. ! ! Output, real Y(N), the product A * X. ! ! Input, real DIAG(NDIAG,IDIAG), the diagonals. ! ! Input, integer NDIAG, the first dimension of array adiag as declared in ! the calling program. ! ! Input, integer IDIAG, the number of diagonals in the matrix. ! ! Input, integer IOFF(IDIAG), the offsets of the diagonals of the matrix: ! diag(i,k) contains the element a(i,i+ioff(k)) of the matrix. ! implicit none integer idiag integer n integer ndiag real ( kind = 8 ) diag(ndiag,idiag) integer i1 integer i2 integer io integer ioff(idiag) integer j integer k real ( kind = 8 ) x(n) real ( kind = 8 ) y(n) y(1:n) = 0.0D+00 do j = 1, idiag io = ioff(j) i1 = max ( 1, 1 - io ) i2 = min ( n, n - io ) do k = i1, i2 y(k) = y(k) + diag(k,j) * x(k+io) end do end do return end subroutine amuxe ( n, x, y, na, ncol, a, ja ) !*****************************************************************************80 ! !! AMUXE multiplies an ELL matrix times a vector. ! ! Discussion: ! ! This routine multiplies a matrix by a vector, where the matrix is stored ! in the ELL Ellpack/Itpack sparse format. ! ! Modified: ! ! 09 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), array of length equal to the column dimension of ! the A matrix. ! ! Input, integer NA, the first dimension of arrays A and JA ! as declared by the calling program. ! ! Input, integer NCOL, the number of active columns in array a. ! (i.e., the number of generalized diagonals in matrix.) ! ! a, ja = the real and integer arrays of the Ellpack/Itpack format ! (a(i,k),k = 1,ncol contains the elements of row i in matrix ! ja(i,k),k = 1,ncol contains their column numbers) ! ! Output, real Y(N), the product A * X. ! implicit none integer n integer na integer ncol real ( kind = 8 ) a(na,ncol) integer i integer j integer ja(na,ncol) real ( kind = 8 ) x(n) real ( kind = 8 ) y(n) y(1:n) = 0.0D+00 do j = 1, ncol do i = 1, n y(i) = y(i) + a(i,j) * x(ja(i,j)) end do end do return end subroutine amuxj ( n, x, y, jdiag, a, ja, ia ) !*****************************************************************************80 ! !! AMUXJ multiplies a JAD matrix times a vector. ! ! Discussion: ! ! This routine multiplies a matrix A times a vector, where A is ! stored in JAD Jagged-Diagonal storage format. ! ! Permutation related to the JAD format is not performed. ! This can be done by: ! ! call permvec ( n, y, y, iperm ) ! ! after the call to AMUXJ. Here IPERM is the permutation produced ! by CSRJAD. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), an array of length equal to the column dimension of ! the A matrix. ! ! Input, integer JDIAG, the number of jagged-diagonals in the data ! structure. ! ! a = real array containing the jagged diagonals of A stored ! in succession (in decreasing lengths) ! ! j = integer array containing the column indices of the ! corresponding elements in a. ! ! ia = integer array containing the lengths of the jagged diagonals ! ! Output, real Y(N), the product A*X. ! implicit none integer n real ( kind = 8 ) a(*) integer ia(*) integer ii integer j integer ja(*) integer jdiag integer k1 integer len real ( kind = 8 ) x(n) real ( kind = 8 ) y(n) y(1:n) = 0.0D+00 do ii = 1, jdiag k1 = ia(ii) - 1 len = ia(ii+1) - k1 - 1 do j = 1, len y(j) = y(j) + a(k1+j) * x(ja(k1+j)) end do 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 NROW, the row dimension of A and B. ! ! Input, integer NCOL, the column dimension of A and B. ! ! Input, integer 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 JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format. ! ! nzmax = integer. 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. 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 real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer ia(nrow+1) integer ib(nrow+1) integer ic(nrow+1) integer ierr integer ii integer iw(ncol) integer ja(*) integer jb(*) integer jc(*) integer jcol integer job integer jpos integer k integer ka integer kb integer len integer 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 NROW, the row dimension of A and B. ! ! Input, integer NCOL, the column dimension of A and B. ! ! Input, integer 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 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. 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. 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 nrow real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer i integer ia(nrow+1) integer ib(nrow+1) integer ic(nrow+1) integer ierr integer j1 integer j2 integer ja(*) integer jb(*) integer jc(*) integer job integer ka integer kamax integer kb integer kbmax integer kc integer ncol integer 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 aplbdg ( nrow, ncol, ja, ia, jb, ib, ndegr, nnz, iw ) !*****************************************************************************80 ! !! APLBDG gets the number of nonzero elements in each row of A + B. ! ! Discussion: ! ! It also reports the total number of nonzero elements in A + B. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of A and B. ! ! Input, integer NCOL, the column dimension of A and B. ! ! 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. ! ! Output, integer NDEGR(NROW), the number of nonzeros in each row ! of the matrix A + B. ! ! Output, integer NNZ, the total number of nonzero elements found ! in A * B. ! ! Workspace, integer IW(NCOL). ! implicit none integer ncol integer nrow integer ia(nrow+1) integer ib(nrow+1) integer ii integer iw(ncol) integer j integer ja(*) integer jb(*) integer jc integer jr integer k integer last integer ldg integer ndegr(nrow) integer nnz iw(1:ncol) = 0 ndegr(1:nrow) = 0 do ii = 1, nrow ldg = 0 ! ! End-of-linked list. ! last = -1 ! ! Row of A. ! do j = ia(ii), ia(ii+1)-1 jr = ja(j) ! ! Add element to the linked list. ! ldg = ldg + 1 iw(jr) = last last = jr end do ! ! Row of B. ! do j = ib(ii), ib(ii+1)-1 jc = jb(j) ! ! Add one element to the linked list. ! if ( iw(jc) == 0 ) then ldg = ldg + 1 iw(jc) = last last = jc end if end do ! ! Done with row II. ! 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 apldia ( nrow, job, a, ja, ia, diag, b, jb, ib, iw ) !*****************************************************************************80 ! !! APLDIA adds a diagonal matrix to a general sparse matrix: B = A + Diag ! ! Discussion: ! ! The column dimension of A is not needed. ! ! The algorithm is in place (b, jb, ib, can be the same as ! a, ja, ia, on entry). See comments for parameter job. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer JOB, job indicator. Job=0 means get array b only ! (i.e. assume that a has already been copied into array b, ! or that algorithm is used in place. ) For all practical ! puposes enter job=0 for an in-place call and job=1 otherwise. ! In case there are missing diagonal elements in A, ! then the option job =0 will be ignored, since the algorithm ! must modify the data structure (i.e. jb, ib) in this ! situation. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, real DIAG(NROW), a diagonal matrix. ! ! on return: ! ! b, ! jb, ! ib = resulting matrix B in compressed sparse row sparse format. ! ! ! iw = integer 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 nrow real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) diag(nrow) integer ia(nrow+1) integer ib(nrow+1) integer icount integer ii integer iw(*) integer j integer ja(*) integer jb(*) integer job integer k integer k1 integer k2 integer ko integer nnz logical test ! ! Copy integer arrays into B's data structure if required. ! if ( job /= 0 ) then nnz = ia(nrow+1)-1 jb(1:nnz) = ja(1:nnz) ib(1:nrow+1) = ia(1:nrow+1) end if ! ! Get positions of diagonal elements in data structure. ! call diapos ( nrow, ja, ia, iw ) ! ! Count number of holes in diagonal and add DIAG elements to ! valid diagonal entries. ! icount = 0 do j = 1, nrow if ( iw(j) == 0 ) then icount = icount + 1 else b(iw(j)) = a(iw(j)) + diag(j) end if end do ! ! If no diagonal elements to insert, return. ! if ( icount == 0 ) then return end if ! ! Shift the nonzero elements if needed, to allow for created ! diagonal elements. ! ko = ib(nrow+1) + icount ! ! Copy rows backward. ! do ii = nrow, 1, -1 ! ! Go through row II. ! k1 = ib(ii) k2 = ib(ii+1) - 1 ib(ii+1) = ko test = ( iw(ii) == 0 ) do k = k2, k1, -1 j = jb(k) if ( test .and. j < ii ) then test = .false. ko = ko - 1 b(ko) = diag(ii) jb(ko) = ii iw(ii) = ko end if ko = ko - 1 b(ko) = a(k) jb(ko) = j end do ! ! The diagonal element has not been added yet. ! if ( test ) then ko = ko - 1 b(ko) = diag(ii) jb(ko) = ii iw(ii) = ko end if end do ib(1) = ko return end subroutine aplsb ( nrow, ncol, a, ja, ia, s, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !*****************************************************************************80 ! !! APLSB performs the matrix linear combination C = A + s * B. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer NCOL, the column dimension of the matrix B. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, real S, scalar factor for B. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format. ! ! nzmax = integer. 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. 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 real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer ia(nrow+1) integer ib(nrow+1) integer ic(nrow+1) integer ierr integer ii integer iw(ncol) integer ja(*) integer jb(*) integer jc(*) integer jcol integer jpos integer k integer ka integer kb integer len integer nzmax real ( kind = 8 ) s 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 c(len) = a(ka) 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 c(len) = s * b(kb) iw(jcol)= len else c(jpos) = c(jpos) + s * b(kb) 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 aplsb1 ( nrow, ncol, a, ja, ia, s, b, jb, ib, c, jc, ic, & nzmax, ierr ) !*****************************************************************************80 ! !! APLSB1 performs the operation C = A + s * B for sorted CSR matrices. ! ! Discussion: ! ! The difference with aplsb is that 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 will not work if any of the two input matrices is not sorted ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of A and B. ! ! Input, integer NCOL, the column dimension of A and B. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format with entries sorted. ! ! Input, real S, a scale factor for B. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format with entries sorted ! ascendly in each row ! ! nzmax = integer. 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. 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 nrow real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer i integer ia(nrow+1) integer ib(nrow+1) integer ic(nrow+1) integer ierr integer j1 integer j2 integer ja(*) integer jb(*) integer jc(*) integer ka integer kamax integer kb integer kbmax integer kc integer ncol integer nzmax real ( kind = 8 ) s 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 c(kc) = a(ka) + s * b(kb) jc(kc) = j1 ka = ka + 1 kb = kb + 1 kc = kc + 1 else if ( j1 < j2 ) then jc(kc) = j1 c(kc) = a(ka) ka = ka + 1 kc = kc + 1 else if ( j2 < j1 ) then jc(kc) = j2 c(kc) = s * b(kb) 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 aplsbt ( nrow, ncol, a, ja, ia, s, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !*****************************************************************************80 ! !! APLSBT performs the matrix sum C = A + B'. ! ! Discussion: ! ! It is important to note that here all of three arrays c, ic, ! and jc are assumed to be of length nnz(c). This is because ! the matrix is internally converted to coordinate "COO" format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of A. This must also be ! the column dimension of B. ! ! Input, integer NCOL, the column dimension of the matrix A. ! This must also be the row dimension of B. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, real S, the scalar factor for B. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format. ! ! nzmax = integer. The length of the arrays c, jc, and ic. ! 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 format. ! ! ierr = integer. serving as error message. ! ierr = 0 means normal return. ! ierr = -1 means that nzmax was < either the number of ! nonzero elements of A or the number of nonzero elements in B. ! 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 real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer ia(nrow+1) integer ib(ncol+1) integer ic(*) integer ierr integer ii integer ipos integer iw(ncol) integer ja(*) integer jb(*) integer jc(*) integer jcol integer jpos integer k integer ka integer len integer ljob integer nnza integer nnzb integer nzmax real ( kind = 8 ) s ierr = 0 iw(1:ncol) = 0 nnza = ia(nrow+1) - 1 nnzb = ib(ncol+1) - 1 len = nnzb if ( nzmax < nnzb .or. nzmax < nnza ) then ierr = -1 return end if ! ! Transpose matrix B into C. ! ljob = 1 ipos = 1 call csrcsc ( ncol, ljob, ipos, b, jb, ib, c, jc, ic ) c(1:len) = c(1:len) * s ! ! The main loop. ! Add rows from 1 through NROW. ! do ii = 1, nrow ! ! IW is used as a system to recognize whether there ! was a nonzero element in C. ! do k = ic(ii), ic(ii+1)-1 iw(jc(k)) = k end do do ka = ia(ii), ia(ii+1)-1 jcol = ja(ka) jpos = iw(jcol) ! ! If fill-in, append in coordinate format to matrix. ! if ( jpos == 0 ) then len = len + 1 if ( nzmax < len ) then ierr = ii return end if jc(len) = jcol ic(len) = ii c(len) = a(ka) else ! ! else, do addition. ! c(jpos) = c(jpos) + a(ka) end if end do do k = ic(ii), ic(ii+1)-1 iw(jc(k)) = 0 end do end do ! ! Convert matrix without fill-ins into coordinate format. ! ljob = 3 call csrcoo ( nrow, ljob, nnzb, c, jc, ic, nnzb, c, ic, jc, ierr ) if ( ierr /= 0 ) then ierr = -ierr end if ! ! Convert the whole thing back to CSR format. ! ljob = 1 call coocsr_inplace ( nrow, len, 1, c, jc, ic, iw ) return end subroutine aplsca ( nrow, a, ja, ia, scal, iw ) !*****************************************************************************80 ! !! APLSCA adds a scalar 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 NROW, the row dimension of the matrix. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, real SCAL, a scalar to be added to the diagonal entries. ! ! on return: ! ! ! a, ! ja, ! ia = matrix A with diagonal elements shifted (or created). ! ! iw = integer 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 nrow real ( kind = 8 ) a(*) integer ia(nrow+1) integer icount integer ii integer iw(*) integer j integer ja(*) integer k integer k1 integer k2 integer ko real ( kind = 8 ) scal 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)) + scal 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) = scal 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) = scal ja(ko) = ii iw(ii) = ko end if end do ia(1) = ko return end subroutine apmbt ( nrow, ncol, job, a, ja, ia, b, jb, ib, c, jc, ic, nzmax, & iw, ierr ) !*****************************************************************************80 ! !! APMBT performs the matrix sum C = A + B' or C = A - B'. ! ! Discussion: ! ! It is important to note that here all of three arrays c, ic, ! and jc are assumed to be of length nnz(c). This is because ! the matrix is internally converted to coordinate "COO" format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of A, which must also be ! the column dimension of B. ! ! Input, integer NCOL, the column dimension of the matrix A, which ! must also be the row dimension of B. ! ! job = integer. if job = -1, apmbt will compute C= A - transp(B) ! (structure + values) ! if job == 1, it will compute C=A+transp(A) ! (structure+ values) ! if job == 0, it will compute the structure of ! C= A+/-transp(B) only (ignoring all real values). ! any other value of job will be treated as job=1 ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! b, ! jb, ! ib = Matrix B in compressed sparse row format. ! ! nzmax = integer. The length of the arrays c, jc, and ic. ! 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 format. ! ! ierr = integer. serving as error message. ! ierr = 0 means normal return. ! ierr = -1 means that nzmax was < either the number of ! nonzero elements of A or the number of nonzero elements in B. ! 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 real ( kind = 8 ) a(*) real ( kind = 8 ) b(*) real ( kind = 8 ) c(*) integer ia(nrow+1) integer ib(ncol+1) integer ic(*) integer ierr integer ii integer ipos integer iw(ncol) integer ja(*) integer jb(*) integer jc(*) integer jcol integer job integer jpos integer k integer ka integer len integer ljob integer nnza integer nnzb integer nzmax logical values values = ( job /= 0 ) ierr = 0 iw(1:ncol) = 0 nnza = ia(nrow+1) - 1 nnzb = ib(ncol+1) - 1 len = nnzb if ( nzmax < nnzb .or. nzmax < nnza ) then ierr = -1 return end if ! ! Transpose matrix B into C. ! ljob = 0 if ( values ) then ljob = 1 end if ipos = 1 call csrcsc ( ncol, ljob, ipos, b, jb, ib, c, jc, ic ) if ( job == -1 ) then c(1:len) = -c(1:len) end if ! ! The main loop. ! do ii = 1, nrow do k = ic(ii), ic(ii+1)-1 iw(jc(k)) = k end do do ka = ia(ii), ia(ii+1)-1 jcol = ja(ka) jpos = iw(jcol) ! ! If fill-in, append in coordinate format to matrix. ! if ( jpos == 0 ) then len = len + 1 if ( nzmax < len ) then ierr = ii return end if jc(len) = jcol ic(len) = ii if ( values ) then c(len) = a(ka) end if else ! ! else do addition. ! if ( values ) then c(jpos) = c(jpos) + a(ka) end if end if end do do k = ic(ii), ic(ii+1)-1 iw(jc(k)) = 0 end do end do ! ! Convert first part of matrix (without fill-ins) into COO format. ! ljob = 2 if ( values ) then ljob = 3 end if call csrcoo ( nrow, ljob, nnzb, c, jc, ic, nnzb ,c, ic, jc, ierr ) if ( ierr /= 0 ) then ierr = -ierr end if ! ! Convert the whole thing back to CSR format. ! ljob = 0 if ( values ) then ljob = 1 end if call coocsr_inplace ( nrow, len, ljob, c, jc, ic, iw ) return end subroutine assmb1 ( u, a, ja, ia, fu, f, node_num, element_num, element_node, & node_code, npe ) !*****************************************************************************80 ! !! ASSMB1 assembles a finite element matrix in the CSR format. ! ! Discussion: ! ! The routine receives all the unassembled local finite element ! matrices and right hand sides, and assembles them into a global ! matrix and right hand side. ! ! The JA and IA arrays are constructed based on the information ! in the element connectivity array ELEMENT_NODE. ! ! Modified: ! ! 01 July 2005 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, real U(ELEMENT_NUM,NPE,NPE), the unassembled local matrices. ! ! Output, real A(*), integer JA(*), IA(NODE_NUM+1), the assembled global ! matrix in CSR (Compressed Sparse Row) format. ! ! Input, real FU(ELEMENT_NUM,NPE), the unassembled right hand sides. ! ! Output, real F(NODE_NUM), the assembled global right hand side. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(NPE,ELEMENT_NUM), the connectivity matrix. ! ELEMENT_NODE(I,J) is the global index of the I-th local node in element J. ! ! Input, integer NODE_CODE(NODE_NUM), boundary information for each node ! with the following meaning: ! * 0, node I is internal; ! * 1, node I is a boundary but not a corner point; ! * 2, node I is a corner point. ! ! Input, integer NPE, the number of nodes per element. ! ! Local parameters: ! ! Workspace, integer IWK(NODE_NUM). ! ! Workspace, integer JWK(NODE_NUM+1). ! ! Integer LOCAL, LOCAL1, LOCAL2, local node numbers. ! ! Integer NODE, NODE1, NODE2, global node numbers. ! implicit none integer element_num integer node_num integer npe real ( kind = 8 ) a(*) integer element integer element_node(npe,element_num) real ( kind = 8 ) f(node_num) real ( kind = 8 ) fu(element_num,npe) integer ia(node_num+1) integer iwk(node_num) integer j integer ja(*) integer jwk(node_num+1) integer k integer k1 integer k2 integer local integer local1 integer local2 integer node integer node_code(node_num) integer node1 integer node2 integer row_last integer row_start real ( kind = 8 ) u(element_num,npe,npe) ! ! Initialize. ! f(1:node_num) = 0.0D+00 ! ! Initialize the pointer arrays. ! ia(1:node_num+1) = 1 jwk(1:node_num+1) = 0 ! ! Count the number of elements (or boundary conditions) where a given ! node occurs. Put this into the IA vector. Then replace this count ! by an incremental count of all the entries preceding it. ! do element = 1, element_num do local = 1, npe node = element_node(local,element) ia(node) = ia(node) + 1 end do end do do node = 1, node_num if ( 1 <= node_code(node) ) then ia(node) = ia(node) + 1 end if end do k1 = ia(1) ia(1) = 1 do node = 2, node_num + 1 k2 = ia(node) ia(node) = ia(node-1) + k1 iwk(node-1) = ia(node-1) - 1 k1 = k2 end do ! ! The assembly loop. ! do element = 1, element_num ! ! The local row number is LOCAL1. ! The global row number is NODE1. ! do local1 = 1, npe node1 = element_node(local1,element) f(node1) = f(node1) + fu(element,local1) ! ! Unpack the row into JWK1. ! row_start = ia(node1) row_last = iwk(node1) do k = row_start, row_last jwk(ja(k)) = k end do ! ! The local column is LOCAL2. ! The global column number is JJ. ! do local2 = 1, npe node2 = element_node(local2,element) k = jwk(node2) if ( k == 0 ) then row_last = row_last + 1 jwk(node2) = row_last ja(row_last) = node2 a(row_last) = u(element,local1,local2) else a(k) = a(k) + u(element,local1,local2) end if end do ! ! Refresh JWK. ! jwk(ja(row_start:row_last)) = 0 iwk(node1) = row_last end do end do return end subroutine assmbo ( nx, nelx, node, ijk, nodcode, x, y, a, ja, ia, f, iwk, & jwk, ierr, xyk ) !*****************************************************************************80 ! !! ASSMBO assembles a finite element matrix. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NX, the number of nodes. ! ! Input, integer NELX, the number of elements. ! ! Input, integer NODE, the number of nodes per element. ! ! Input, integer IJK(NODE,NELX), lists the nodes that make up ! each element. ! ! Output, real A(*), integer JA(*), IA(NX+1), the assembled matrix in CSR ! Compressed Sparse Row format. ! ! f = right hand side (global load vector) ! ! nodcode= boundary information list for each node with the ! following meaning: ! nodcode(i) = 0 --> node i is internal ! nodcode(i) = 1 --> node i is a boundary but not a corner point ! nodcode(i) = 2 --> node i is a corner point (corner points ! ! mp = group material number for each element. ! ! x,y = real arrays containing the x and y coordinates ! resp. of the nodes. ! ! iwk,jwk = two integer work arrays. ! ierr = error message integer . ! ierr = 0 --> normal return ! ierr = 1 --> negative area encountered (due to bad ! numbering of nodes of an element- see ! message printed in unit iout ). ! iout = output unit (not used here). ! ! xyk = routine defining the material properties at each ! element. Form: ! call xyk(nel,xyke,x,y,ijk,node) with on return ! xyke = material constant matrices. ! for each element nel, xyke(1,nel),xyke(2,nel) ! and xyke(3,nel) represent the constants ! K11, K22, and K12 at that element. ! implicit none integer node integer nx real ( kind = 8 ) a(*) real ( kind = 8 ) det real ( kind = 8 ) f(*) real ( kind = 8 ) fe(3) integer i integer ia(nx+1) integer ierr integer ii integer ijk(node,*) integer ilast integer irowst integer iwk(*) integer j integer ja(*) integer jj integer jwk(*) integer k integer ka integer kb integer knod integer ksav integer ksavn integer nel integer nelx integer nodcode(*) real ( kind = 8 ) ske(3,3) real ( kind = 8 ) x(*) real ( kind = 8 ) xe(3) external xyk real ( kind = 8 ) xyke(2,2) real ( kind = 8 ) y(*) real ( kind = 8 ) ye(3) ! ! Maximum number of nonzeros allowed = 200 ! ! Initialize. ! f(1:nx) = 0.0D+00 ! ! Initialize the pointer arrays. ! ia(1:nx+1) = 1 jwk(1:nx+1) = 0 do k = 1, nelx do j = 1, node knod = ijk(j,k) ia(knod) = ia(knod) + 1 end do end do do k = 1, nx if ( 1 <= nodcode(k) ) then ia(k) = ia(k) + 1 end if end do ksav = ia(1) ia(1) = 1 do j = 2, nx+1 ksavn = ia(j) ia(j) = ia(j-1) + ksav iwk(j-1) = ia(j-1)-1 ksav = ksavn end do ! ! The main loop. ! do nel = 1, nelx ! ! Get coordinates of nodal points. ! do i = 1, node j = ijk(i,nel) xe(i) = x(j) ye(i) = y(j) end do ! ! Compute determinant. ! det = xe(2) * ( ye(3) - ye(1) ) & + xe(3) * ( ye(1) - ye(2) ) & + xe(1) * ( ye(2) - ye(3) ) ! ! Set material properties. ! call xyk ( nel, xyke, x, y, ijk, node ) ! ! Construct element stiffness matrix. ! ierr = 0 call estif3 ( nel, ske, fe, det, xe, ye, xyke, ierr ) if ( ierr /= 0 ) then return end if ! ! Assemble: add element stiffness matrix to global matrix. ! do ka = 1, node ii = ijk(ka,nel) f(ii) = f(ii) + fe(ka) ! ! Unpack row into JWK1. ! irowst = ia(ii) ilast = iwk(ii) do k = irowst, ilast jwk(ja(k)) = k end do do kb = 1, node ! ! Column number = JJ. ! jj = ijk(kb,nel) k = jwk(jj) if ( k == 0 ) then ilast = ilast + 1 jwk(jj) = ilast ja(ilast) = jj a(ilast) = ske(ka,kb) else a(k) = a(k) + ske(ka,kb) end if end do ! ! Refresh JWK. ! do k = irowst, ilast jwk(ja(k)) = 0 end do iwk(ii) = ilast end do end do return end subroutine atmux ( n, x, y, a, ja, ia ) !*****************************************************************************80 ! !! ATMUX computes A' * x for a CSR matrix A. ! ! Discussion: ! ! This routine multiplies the transpose of a matrix by a vector when the ! original matrix is stored in compressed sparse row storage. Can also be ! viewed as the product of a matrix by a vector when the original ! matrix is stored in the compressed sparse column format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real X(*), an array whose length is equal to the ! column dimension of A. ! ! Output, real Y(N), the product A' * X. ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! implicit none integer n real ( kind = 8 ) a(*) integer i integer ia(*) integer ja(*) integer k real ( kind = 8 ) x(*) real ( kind = 8 ) y(n) y(1:n) = 0.0D+00 do i = 1, n do k = ia(i), ia(i+1)-1 y(ja(k)) = y(ja(k)) + x(i) * a(k) end do end do return end subroutine blkchk ( nrow, ja, ia, nblk, imsg ) !*****************************************************************************80 ! !! BLKCHK checks whether the input matrix is a block matrix. ! ! Discussion: ! ! This routine checks whether the input matrix is a block matrix with block ! size of NBLK. A block matrix is one which is comprised of small square ! dense blocks. If there are zero elements within the square blocks and the ! data structure takes them into account then blkchk may fail to find the ! correct block size. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, JA(*), IA(NROW+1), the matrix information (but no values) in CSR ! Compressed Sparse Row format. ! ! Input, integer NBLK, the block value to be checked. ! ! Output, integer IMSG, a message with the following meaning: ! 0 : the output value of NBLK is a correct block size. ! -1 : NBLK does not divide NROW; ! -2 : a starting element in a row is at wrong position ! (j /= mult*nblk +1 ); ! -3 : NBLK does not divide a row length; ! -4 : an element is isolated outside a block or two rows in same ! group have different lengths ! implicit none integer nrow integer i integer i1 integer ia(nrow+1) integer ii integer imsg integer irow integer j integer j2 integer ja(*) integer jstart integer k integer len integer lena integer nblk integer nr ! ! First part of code will find candidate block sizes. ! This is not guaranteed to work, so a check is done at the end. ! The criterion used here is a simple one: ! scan rows and determine groups of rows that have the same length ! and such that the first column number and the last column number ! are identical. ! imsg = 0 if ( nblk <= 1 ) then return end if nr = nrow / nblk if ( nr * nblk /= nrow ) then imsg = -1 return end if ! ! The main loop. ! irow = 1 do ii = 1, nr ! ! I1 = starting position for group of NBLK rows in original matrix. ! i1 = ia(irow) j2 = i1 ! ! LENA = length of each row in that group in the original matrix. ! lena = ia(irow+1) - i1 ! ! LEN = length of each block-row in that group in the output matrix. ! len = lena / nblk if ( len * nblk /= lena ) then imsg = -3 return end if ! ! For each row. ! do i = 1, nblk irow = irow + 1 if ( ia(irow) - ia(irow-1) /= lena ) then imsg = -4 return end if ! ! For each block. ! do k = 0, len-1 jstart = ja(i1+nblk*k) - 1 if ( ( jstart / nblk ) * nblk /= jstart ) then imsg = -2 return end if ! ! For each column. ! do j = 1, nblk if ( jstart + j /= ja(j2) ) then imsg = -4 end if j2 = j2 + 1 end do end do end do end do return end subroutine blkfnd ( nrow, ja, ia, nblk ) !*****************************************************************************80 ! !! BLKFND determines whether the input matrix has a block structure. ! ! Discussion: ! ! If the matrix has a block structure, this routine finds the block ! size. A block matrix is one which is comprised of small square ! dense blocks. If there are zero elements within the square blocks ! and the original data structure takes these zeros into account ! then this routine may fail to find the correct block size. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, real JA(*), IA(NROW+1), the matrix information (but not the ! values) in CSR Compressed Sparse Row format. ! ! Output, integer NBLK, the block value that was found. ! implicit none integer nrow integer i integer i1 integer i2 integer ia(nrow+1) integer iblk integer imsg integer irow integer ja(*) integer jf integer jfirst integer jl integer jlast integer jrow integer len integer len0 integer minlen integer nblk ! ! The first part of code will find candidate block sizes. ! The criterion used here is a simple one: scan rows and determine groups ! of rows that have the same length and such that the first column ! number and the last column number are identical. ! minlen = ia(2) - ia(1) irow = 1 do i = 2, nrow len = ia(i+1) - ia(i) if ( len < minlen ) then minlen = len irow = i end if end do ! ! Candidates are all dividers of MINLEN. ! nblk = 1 if ( minlen <= 1 ) then return end if do iblk = minlen, 1, -1 if ( mod ( minlen, iblk ) /= 0 ) then cycle end if len = ia(2) - ia(1) len0 = len jfirst = ja(1) jlast = ja(ia(2)-1) do jrow = irow+1, irow+nblk-1 i1 = ia(jrow) i2 = ia(jrow+1) - 1 len = i2 + 1 - i1 jf = ja(i1) jl = ja(i2) if ( len /= len0 .or. jf /= jfirst .or. jl /= jlast ) then go to 99 end if end do ! ! Check for this candidate. ! call blkchk ( nrow, ja, ia, iblk, imsg ) ! ! Block size found. ! if ( imsg == 0 ) then nblk = iblk return end if 99 continue end do return end subroutine bndcsr ( n, abd, nabd, lowd, ml, mu, a, ja, ia, len, ierr ) !*****************************************************************************80 ! !! BNDCSR converts Banded Linpack format to Compressed Sparse Row format. ! ! Discussion: ! ! The matrix values found to be equal to zero ! (actual test: if (abd(...) == 0.0) are removed. ! ! The resulting may not be identical to a CSR matrix ! originally transformed to a BND format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, integer NABD, the first dimension of ABD. ! ! abd = real array containing the values of the matrix stored in ! banded form. The j-th column of abd contains the elements ! of the j-th column of the original matrix,comprised in the ! band ( i in (j-ml,j+mu) ) with the lowest diagonal located ! in row lowd (see below). ! ! lowd = integer. this should be set to the row number in abd where ! the lowest diagonal (leftmost) of A is located. ! lowd should be s.t. ( 1 <= lowd <= nabd). ! The routines dgbco, ... of linpack use lowd=2*ml+mu+1. ! ! ml = integer. equal to the bandwidth of the strict lower part of A ! mu = integer. equal to the bandwidth of the strict upper part of A ! thus the total bandwidth of A is ml+mu+1. ! if ml+mu+1 is found to be larger than nabd then an error ! message is set. see ierr. ! ! len = integer. length of arrays a and ja. bndcsr will stop if the ! length of the arrays a and ja is insufficient to store the ! matrix. see ierr. ! ! on return: ! ! Output, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! ! lowd = if on entry lowd was zero then lowd is reset to the default ! value ml+mu+l. ! ! ierr = integer. used for error message output. ! ierr == 0 :means normal return ! ierr == -1 : means invalid value for lowd. ! ierr > 0 : means that there was not enough storage in a and ja ! for storing the ourput matrix. The process ran out of space ! (as indicated by len) while trying to fill row number ierr. ! This should give an idea of much more storage might be required. ! Moreover, the first irow-1 rows are correctly filled. ! implicit none integer n integer nabd real ( kind = 8 ) a(*) real ( kind = 8 ) abd(nabd,*) integer i integer ia(n+1) integer ierr integer irow integer j integer ja(*) integer ko integer len integer lowd integer ml integer mu real ( kind = 8 ) t ierr = 0 if ( nabd < lowd .or. lowd <= 0 ) then ierr = -1 return end if ko = 1 ia(1) = 1 do irow = 1, n i = lowd do j = irow-ml, irow+mu if ( j <= 0 ) then go to 19 end if if ( n < j ) then go to 21 end if t = abd(i,j) if ( t /= 0.0D+00 ) then if ( len < ko ) then ierr = irow return end if a(ko) = t ja(ko) = j ko = ko + 1 end if 19 continue i = i - 1 end do 21 continue ia(irow+1) = ko end do return end subroutine bound ( nx, nelx, ijk, nodcode, node, n_int, iperm, & x, y, wk, iwk ) !*****************************************************************************80 ! !! BOUND counts the number of boundary points. ! ! Discussion: ! ! It also reorders the points in such a way that the boundary nodes ! are last. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NX, the number of nodes. ! ! Input, integer NELX, the number of elements. ! ! Input/output, integer IJK(NODE,NELX), lists the nodes that make up ! each element. On output, IJK has been updated. ! ! nodcode, node: see other routines ! ! Output, integer N_INT, the number of points on the boundary. ! ! iperm = permutation array from old ordering to new ordering, ! ! iwk = reverse permutation array or return. ! wk = real work array ! On return ! x, y, nodecode, are permuted ! ijk is updated according to new oerdering. ! n_int = number of interior points. ! implicit none integer node integer ijk(node,*) integer iperm(*) integer iwk(*) integer j integer k integer knod integer n_int integer nbound integer nel integer nelx integer nodcode(*) integer nx real ( kind = 8 ) wk(*) real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) ! ! Maximum number of nonzeros allowed = 200 ! ! Put all boundary points at the end, backwards. ! n_int = 1 nbound = nx do j = 1, nx if ( nodcode(j) == 0 ) then iperm(n_int) = j n_int = n_int + 1 else iperm(nbound) = j nbound = nbound - 1 end if end do n_int = n_int - 1 ! ! Permute X's. ! wk(1:nx) = x(1:nx) x(1:nx) = wk(iperm(1:nx)) ! ! Permute the Y's. ! wk(1:nx) = y(1:nx) y(1:nx) = wk(iperm(1:nx)) ! ! Permute the boundary information. ! iwk(1:nx) = nodcode(1:nx) do k = 1, nx nodcode(k) = iwk(iperm(k)) end do ! ! Get the reverse permutation. ! do k = 1, nx iwk(iperm(k)) = k end do ! ! Update the element connectivity matrix. ! do nel = 1, nelx do j = 1, node knod = ijk(j,nel) ijk(j,nel) = iwk(knod) end do end do return end subroutine bsort2 ( w, ind, n, ncut ) !*****************************************************************************80 ! !! BSORT2 returns the NCUT largest elements of an array, using bubble sort. ! ! Discussion: ! ! This routine carries out a simple bubble sort for getting the NCUT largest ! elements in modulus, in array W. IND is sorted accordingly. ! (Ought to be replaced by a more efficient sort especially ! if NCUT is not that small). ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! implicit none integer n integer i integer ind(*) integer iswp integer j integer ncut logical test real ( kind = 8 ) w(n) real ( kind = 8 ) wswp i = 1 do test = .false. do j = n-1, i, -1 if ( abs ( w(j) ) < abs ( w(j+1) ) ) then ! ! Swap. ! wswp = w(j) w(j) = w(j+1) w(j+1) = wswp ! ! Reorder the original ind array accordingly. ! iswp = ind(j) ind(j) = ind(j+1) ind(j+1) = iswp ! ! Set indicator that sequence is still unsorted. ! test = .true. end if end do i = i + 1 if ( .not. test .or. ncut < i ) then exit end if end do return end subroutine bsrcsr ( n, nblk, na, a, ja, ia, ao, jao, iao ) !*****************************************************************************80 ! !! BSRCSR converts Block Sparse Row to Compressed Sparse Row (CSR) format. ! ! Discussion: ! ! This routine converts a matrix stored in block-reduced ! a, ja, ia format to the general sparse row a, ja, ia format. ! A matrix that has a block structure is a matrix whose entries ! are blocks of the same size nblk (e.g. 3 x 3). Then it is often ! preferred to work with the reduced graph of the matrix, i.e., ! Instead of storing one element at a time one can store the whole ! block. In this storage scheme a row of the array a will ! hold the nblk**2 entries of a block. ! ! This code is not "in place". ! ! general picture: (nblk = 2) ! --- A --- --- JA -- -- IA -- ! A= x x x x 1st block in block row 1 x x ! x x x x 2-nd block in block row 1 x ! . . . . . ! x x x x last block in block row 1 x ! ------- --- ! x x x x 1st block in block row 2 x x ! x x x x 2-nd block in block row 2 x ! . . . . x ! x x x x last block in block row 2 x ! ------- --- ! ....... ... . ! ------- --- ! x x x x 1st block in block row n/nblk x x ! x x x x 2-nd block in block row n/nblk x ! . . . . x ! x x x x last block in block row n/nblk x ! ------- --- ! end + 1 x ! ! ! example with nblk = 2: ! ! ! 1 2 0 0 3 4 ! 5 6 0 0 7 8 ! 0 0 9 10 11 12 ! 0 0 13 14 15 16 ! 17 18 0 0 0 0 ! 22 23 0 0 0 0 ! THEN: ! ! ---- A ---- -- JA -- -- IA -- !- ! 1 5 2 6 Block row 1 (2 block matrices) | 1 <--- | 1 ! 3 7 4 8 | 5 | ! ------------ | -- | ! 9 13 10 14 block row 2 (2 block matrices) | 3 <--- | 3 ! 11 15 12 16 | 5 | ! ------------ | -- | ! 17 22 18 23 Block row 3 (1 block matrix) | 1 <--- | 5 ! ------------ | -- | ! end+1 <--- | 6 ! ! JA = 1 5 | 3 5 | 1 column numbers of (1,1) entries of blocks ! IA = 1 3 5 6 pointers to beginnings of BLOCK-rows ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! nblk = integer equal to the dimension of each block. ! nblk must divide n. ! ! na = first dimension of array a as declared in calling program ! ! a = real array containing the values of the matrix. For details ! on the format see below. Each row of a contains the nblk x nblk ! block matrix unpacked column-wise (this allows the user to ! declare the array a as a(na,nblk,nblk) on entry if desired). ! the block rows are stored in sequence just as for the compressed ! sparse row format. ! ! ja = integer array of length n/nblk. ja(k) contains the column index ! of the leading element, i.e., the element (1,1) of the block ! that is held in the row a(k,*) of the value array. ! ! ia = integer array of length n/nblk+1. ia(i) points to the beginning ! of block row number i in the arrays a and ja. ! ! Output, real AO(*), JAO(*), IAO(N+1), the matrix in CSR ! Compressed Sparse Row format. ! implicit none integer n integer na real ( kind = 8 ) a(na,*) real ( kind = 8 ) ao(*) integer i integer i1 integer i2 integer ia(*) integer iao(n+1) integer ii integer ij integer irow integer j integer ja(*) integer jao(*) integer jj integer jst integer k integer krow integer nblk integer nr ! ! Get the IA, JA data structure for output matrix ! nr = n / nblk iao(1:n+1) = 0 irow = 0 krow = 1 do ii = 1, nr ! ! NR is the dimension of the reduced matrix. ! i1 = ia(ii) i2 = ia(ii+1) - 1 ! ! Create NBLK rows for each K. ! do i = 1, nblk do k = i1, i2 jst = ja(k) - 1 do j = 1, nblk ij = ( j - 1 ) * nblk + i ao(krow) = a(k,ij) jao(krow) = jst + j krow = krow + 1 end do end do iao(irow+i) = krow end do irow = irow + nblk end do do jj = 1, n j = n - jj + 1 iao(j+1) = iao(j) end do iao(1) = 1 return end subroutine bsten ( nx, ny, nz, kx, ky, kz, nfree, stencil, h ) !*****************************************************************************80 ! !! BSTEN calculates block stencil values. ! ! Discussion: ! ! This routine calculates the correct block-stencil values for ! a centered difference discretization of the elliptic operator ! (block version of stencil) ! ! L u = delx( a delx u ) + dely ( b dely u) + delz ( c delz u ) + ! d delx ( u ) + e dely (u) + f delz( u ) + g u ! ! For 2-D problems the discretization formula that is used is: ! ! h**2 * Lu == a(i+1/2,j) * {u(i+1,j) - u(i,j)} + ! a(i-1/2,j) * {u(i-1,j) - u(i,j)} + ! b(i,j+1/2) * {u(i,j+1) - u(i,j)} + ! b(i,j-1/2) * {u(i,j-1) - u(i,j)} + ! (h/2) * d(i,j) * {u(i+1,j) - u(i-1,j)} + ! (h/2) * e(i,j) * {u(i,j+1) - u(i,j-1)} + ! (h/2) * e(i,j) * {u(i,j+1) - u(i,j-1)} + ! (h**2) * g(i,j) * u(i,j) ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! implicit none real ( kind = 8 ) cntr(225) real ( kind = 8 ) coeff(225) real ( kind = 8 ) h real ( kind = 8 ) h2 real ( kind = 8 ) hhalf integer k integer kx integer ky integer kz integer nfree integer nfree2 integer nx integer ny integer nz real ( kind = 8 ) stencil(7,*) real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) z if ( 15 < nfree ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BSTEN - FATAL ERROR' write ( *, '(a)' ) ' Input value of NFREE is greater than 15.' stop end if nfree2 = nfree * nfree cntr(1:nfree2) = 0.0D+00 stencil(1:7,1:nfree2) = 0.0D+00 hhalf = h * 0.5D+00 h2 = h * h x = h * real ( kx, kind = 8 ) y = h * real ( ky, kind = 8 ) z = h * real ( kz, kind = 8 ) ! ! Differentiation with respect to X: ! call afunbl ( nfree, x + hhalf, y, z, coeff ) do k = 1, nfree2 stencil(3,k) = stencil(3,k) + coeff(k) cntr(k) = cntr(k) + coeff(k) end do call afunbl ( nfree, x - hhalf, y, z, coeff ) do k = 1, nfree2 stencil(2,k) = stencil(2,k) + coeff(k) cntr(k) = cntr(k) + coeff(k) end do call dfunbl ( nfree, x, y, z, coeff ) do k = 1, nfree2 stencil(3,k) = stencil(3,k) + coeff(k) * hhalf stencil(2,k) = stencil(2,k) - coeff(k) * hhalf end do if ( ny <= 1 ) then go to 99 end if ! ! Differentiation with respect to Y: ! call bfunbl ( nfree, x, y + hhalf, z, coeff ) do k = 1, nfree2 stencil(5,k) = stencil(5,k) + coeff(k) cntr(k) = cntr(k) + coeff(k) end do call bfunbl ( nfree, x, y - hhalf, z, coeff ) do k = 1, nfree2 stencil(4,k) = stencil(4,k) + coeff(k) cntr(k) = cntr(k) + coeff(k) end do call efunbl ( nfree, x, y, z, coeff ) do k = 1, nfree2 stencil(5,k) = stencil(5,k) + coeff(k) * hhalf stencil(4,k) = stencil(4,k) - coeff(k) * hhalf end do ! ! Differentiation with respect to Z: ! if ( 1 < nz ) then call cfunbl ( nfree, x, y, z + hhalf, coeff ) do k = 1, nfree2 stencil(7,k) = stencil(7,k) + coeff(k) cntr(k) = cntr(k) + coeff(k) end do call cfunbl ( nfree, x, y, z - hhalf, coeff ) do k = 1, nfree2 stencil(6,k) = stencil(6,k) + coeff(k) cntr(k) = cntr(k) + coeff(k) end do call ffunbl ( nfree, x, y, z, coeff ) do k = 1, nfree2 stencil(7,k) = stencil(7,k) + coeff(k) * hhalf stencil(6,k) = stencil(6,k) - coeff(k) * hhalf end do end if ! ! Discretization of product by G: ! 99 continue call gfunbl ( nfree, x, y, z, coeff ) do k = 1, nfree2 stencil(1,k) = h2*coeff(k) - cntr(k) end do return end subroutine checkref ( nx, nelx, ijk, node, nodcode, nbound, nxnew, nelxnew ) !*****************************************************************************80 ! !! CHECKREF returns the expected number of new nodes and elements. ! ! Discussion: ! ! These numbers indicate the number of new nodes and elements that may ! be expected if routine REFALL is applied once to the current grid. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NX, the number of nodes. ! ! Input, integer NELX, the number of elements. ! ! Input, integer IJK(NODE,NELX), lists the nodes that make up ! each element. ! ! nbound = number of boundary points on entry - enter zero if ! unknown ! ! nodcode= boundary information list for each node with the ! following meaning: ! nodcode(i) = 0 --> node i is internal ! nodcode(i) = 1 --> node i is a boundary but not a corner point ! nodcode(i) = 2 --> node i is a corner point. ! ! nxnew = new number of nodes if refall were to be applied ! nelxnew = same for nelx. ! implicit none integer node integer nx integer ijk(node,*) integer j integer nbound integer nelx integer nelxnew integer nodcode(nx) integer nxnew nelxnew = nelx * 4 ! ! Count the boundary nodes. ! if ( nbound == 0 ) then do j = 1, nx if ( 1 <= nodcode(j) ) then nbound = nbound + 1 end if end do end if ! ! Number of edges = ( 3 * ( number of elements ) + number of bound nodes ) / 2 ! nxnew = nx + ( 3 * nelx + nbound ) / 2 nbound = 2 * nbound return end subroutine chkelmt ( nx, x, y, nelx, ijk, node ) !*****************************************************************************80 ! !! CHKELMT checks the labeling within each element and reorders if necessary. ! ! Discussion: ! ! If the nodes are not correctly ordered, this routine reorders them. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NX, the number of nodes. ! ! Input, real X(*), Y(*), the coordinates of the nodes. ! ! Input, integer NELX, the number of elements. ! ! Input, integer IJK(NODE,NELX), lists the nodes that make up ! each element. ! ! Input, integer NODE, the number of nodes per element. ! implicit none integer node real ( kind = 8 ) det integer ijk(node,*) integer j integer nel integer nelx integer nx real ( kind = 8 ) x(*) real ( kind = 8 ) y(*) do nel = 1, nelx det = x(ijk(2,nel)) * ( y(ijk(3,nel)) - y(ijk(1,nel)) ) + & x(ijk(3,nel)) * ( y(ijk(1,nel)) - y(ijk(2,nel)) ) + & x(ijk(1,nel)) * ( y(ijk(2,nel)) - y(ijk(3,nel)) ) ! ! If the determinant is negative, switch the last two nodes of the element. ! if ( det < 0.0D+00 ) then j = ijk(2,nel) ijk(2,nel) = ijk(3,nel) ijk(3,nel) = j end if end do return end subroutine cnrms ( nrow, nrm, a, ja, ia, diag ) !*****************************************************************************80 ! !! CNRMS gets the norms of each column of A. ! ! Discussion: ! ! There is a choice of three norms. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer NRM, choosed the norm: ! 1, means 1-norm, ! 2, means the 2-nrm, ! 0, means max norm ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real ( kind = 8 ) DIAG(NROW), the row norms. ! implicit none integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) diag(nrow) integer ia(nrow+1) integer ii integer j integer ja(*) integer k integer k1 integer k2 integer nrm diag(1:nrow) = 0.0D+00 do ii = 1, nrow k1 = ia(ii) k2 = ia(ii+1) - 1 do k = k1, k2 j = ja(k) ! ! Update the norm of each column. ! if ( nrm == 0 ) then diag(j) = max ( diag(j), abs ( a(k) ) ) else if ( nrm == 1 ) then diag(j) = diag(j) + abs ( a(k) ) else diag(j) = diag(j) + a(k)**2 end if end do end do if ( nrm /= 2 ) then return end if do k = 1, nrow diag(k) = sqrt ( diag(k) ) end do return end subroutine coocsr_inplace ( n, nnz, job, a, ja, ia, iwk ) !*****************************************************************************80 ! !! COOCSR_INPLACE converts COO to CSR in place. ! ! Discussion: ! ! This routine converts a matrix stored in coordinate format into ! the CSR format. The conversion is done in place in that the arrays ! a,ja,ia of the result are overwritten onto the original arrays. ! ! The entries of the output matrix are not sorted (the column ! indices in each are not in increasing order) use COOCSR ! if you want them sorted. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, integer NNZ, the number of nonzero elements in A. ! ! Input, integer JOB. When JOB = 1, the real values in A are ! filled. Otherwise A is not touched and the structure of the ! array only (i.e. JA, IA) is obtained. ! ! Input/output, real A(NNZ). On input, the matrix numeric values, ! stored in the COO format. On output, the numeric values, stored ! in CSR format. ! ! ja = integer array of length nnz containing the column positions ! of the corresponding elements in a. ! ! ia = integer array of length nnz containing the row positions ! of the corresponding elements in a. ! ! iwk = integer work array of length n. ! ! on return: ! ! Output, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! implicit none integer n integer nnz real ( kind = 8 ) a(*) integer i integer ia(nnz) integer inext integer init integer ipos integer iwk(n) integer j integer ja(nnz) integer jnext integer job integer k real ( kind = 8 ) t real ( kind = 8 ) tnext logical values values = (job == 1) ! ! Find pointer array for resulting matrix. ! iwk(1:n+1) = 0 do k = 1, nnz i = ia(k) iwk(i+1) = iwk(i+1) + 1 end do iwk(1) = 1 do i = 2, n iwk(i) = iwk(i-1) + iwk(i) end do ! ! Loop for a cycle in chasing process. ! init = 1 k = 0 5 continue if ( values ) then t = a(init) end if i = ia(init) j = ja(init) ia(init) = -1 6 continue k = k + 1 ! ! Current row number is I. Determine where to go. ! ipos = iwk(i) ! ! Save the chased element. ! if ( values ) then tnext = a(ipos) end if inext = ia(ipos) jnext = ja(ipos) ! ! Then occupy its location. ! if ( values ) then a(ipos) = t end if ja(ipos) = j ! ! Update pointer information for next element to come in row I. ! iwk(i) = ipos + 1 ! ! Determine the next element to be chased. ! if ( ia(ipos) < 0 ) then go to 65 end if t = tnext i = inext j = jnext ia(ipos) = -1 if ( k < nnz ) then go to 6 end if go to 70 65 continue init = init + 1 if ( nnz < init ) then go to 70 end if if ( ia(init) < 0 ) then go to 65 end if ! ! Restart chasing. ! go to 5 70 continue ia(1) = 1 ia(2:n+1) = iwk(1:n) return end subroutine coocsr ( nrow, nnz, a, ir, jc, ao, jao, iao ) !*****************************************************************************80 ! !! COOCSR converts COO to CSR. ! ! Discussion: ! ! This routine converts a matrix that is stored in COO coordinate format ! a, ir, jc into a CSR row general sparse ao, jao, iao format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, integer NNZ, the number of nonzero elements in the matrix. ! ! a, ! ir, ! jc = matrix in coordinate format. a(k), ir(k), jc(k) store the nnz ! nonzero elements of the matrix with a(k) = actual real value of ! the elements, ir(k) = its row number and jc(k) = its column ! number. The order of the elements is arbitrary. ! ! on return: ! ! ir is destroyed ! ! Output, real AO(*), JAO(*), IAO(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! implicit none integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) ao(*) integer i integer iad integer iao(nrow+1) integer ir(*) integer j integer jao(*) integer jc(*) integer k integer k0 integer nnz real ( kind = 8 ) x iao(1:nrow+1) = 0 ! ! Determine the row lengths. ! do k = 1, nnz iao(ir(k)) = iao(ir(k)) + 1 end do ! ! The starting position of each row. ! k = 1 do j = 1, nrow+1 k0 = iao(j) iao(j) = k k = k + k0 end do ! ! Go through the structure once more. Fill in output matrix. ! do k = 1, nnz i = ir(k) j = jc(k) x = a(k) iad = iao(i) ao(iad) = x jao(iad) = j iao(i) = iad + 1 end do ! ! Shift back IAO. ! do j = nrow, 1, -1 iao(j+1) = iao(j) end do iao(1) = 1 return end subroutine cooell ( n, nnz, a, ja, ia, ac, jac, nac, ner, ncmax, ierr ) !*****************************************************************************80 ! !! COOELL converts coordinate format to Ellpack/Itpack format. ! ! Discussion: ! ! This routine takes a sparse matrix in coordinate format and ! converts it into the Ellpack/Itpack storage. ! ! Example: ! ! ( 11 0 13 0 0 0 ) ! | 21 22 0 24 0 0 | ! | 0 32 33 0 35 0 | ! A = | 0 0 43 44 0 46 | ! | 51 0 0 54 55 0 | ! ( 61 62 0 0 65 66 ) ! ! Coordinate storage scheme: ! ! A = (11,22,33,44,55,66,13,21,24,32,35,43,46,51,54,61,62,65) ! IA = (1, 2, 3, 4, 5, 6, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 6 ) ! JA = ( 1, 2, 3, 4, 5, 6, 3, 1, 4, 2, 5, 3, 6, 1, 4, 1, 2, 5) ! ! Ellpack/Itpack storage scheme: ! ! ( 11 13 0 0 ) ( 1 3 * * ) ! | 22 21 24 0 | | 2 1 4 * | ! AC = | 33 32 35 0 | JAC = | 3 2 5 * | ! | 44 43 46 0 | | 4 3 6 * | ! | 55 51 54 0 | | 5 1 4 * | ! ( 66 61 62 65 ) ( 6 1 2 5 ) ! ! Note: * means that you can store values from 1 to 6 (1 to n, where ! n is the order of the matrix) in that position in the array. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Ernest E. Rothman, Cornell Theory Center ! ! Reference: ! ! David Kincaid, T C Oppe, J R Respess, D M Young, ! ITPACKV 2C User's Guide, ! Technical Report CNA-191. ! Center for Numerical Analysis, ! University of Texas at Austin, 1984. ! ! Engineering and Scientific Subroutine Library; ! Guide and Reference; ! Release 3 (SC23-0184-3), pages 79-86. ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input, integer NNZ, the number of nonzero elements in the sparse matrix. ! ! Input, integer NCA, the first dimension of output arrays CA and JAC. ! ! A(NNZ) - Real array. ! Stored entries of the sparse matrix A. ! NNZ is the number of nonzeros. ! ! IA(NNZ) - Integer array. ! Pointers to specify rows for the stored nonzero entries ! in A. ! ! JA(NNZ) - Integer array. ! Pointers to specify columns for the stored nonzero ! entries in A. ! ! NER - Integer. Must be set greater than or equal to the maximum ! number of nonzeros in any row of the sparse matrix. ! ! OUTPUT PARAMETERS ! ! AC(NAC,*) - Real array. ! Stored entries of the sparse matrix A in compressed ! storage mode. ! ! JAC(NAC,*) - Integer array. ! Contains the column numbers of the sparse matrix ! elements stored in the corresponding positions in ! array AC. ! ! NCMAX - Integer. Equals the maximum number of nonzeros in any ! row of the sparse matrix. ! ! IERR - Error parameter is returned as zero on successful ! execution of the subroutin= lowd ! ! example [from linpack ]: if the original matrix is ! ! 11 12 13 0 0 0 ! 21 22 23 24 0 0 ! 0 32 33 34 35 0 original banded matrix ! 0 0 43 44 45 46 ! 0 0 0 54 55 56 ! 0 0 0 0 65 66 ! ! then n = 6, ml = 1, mu = 2. lowd should be >= 4 (=ml+mu+1) and ! if lowd = 5 for example, abd should be: ! ! untouched --> x x x x x x ! * * 13 24 35 46 ! * 12 23 34 45 56 resulting abd matrix in banded ! 11 22 33 44 55 66 format ! row lowd--> 21 32 43 54 65 * ! ! * = not used ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! job = integer. if job=1 then the values of the lower bandwith ml ! and the upper bandwidth mu are determined internally. ! otherwise it is assumed that the values of ml and mu ! are the correct bandwidths on input. See ml and mu below. ! ! nabd = integer. first dimension of array abd. ! ! lowd = integer. this should be set to the row number in abd where ! the lowest diagonal (leftmost) of A is located. ! lowd should be ( 1 <= lowd <= nabd). ! if it is not known in advance what lowd should be ! enter lowd = 0 and the default value lowd = ml+mu+1 ! will be chosen. Alternative: call routine getbwd from unary ! first to detrermione ml and mu then define lowd accordingly. ! (Note: the banded solvers in linpack use lowd=2*ml+mu+1. ) ! ! ml = integer. equal to the bandwidth of the strict lower part of A ! mu = integer. equal to the bandwidth of the strict upper part of A ! thus the total bandwidth of A is ml+mu+1. ! if ml+mu+1 is found to be larger than lowd then an error ! flag is raised (unless lowd = 0). see ierr. ! ! note: ml and mu are assumed to have the correct bandwidth values ! as defined above if job is set to zero on entry. ! ! on return: ! ! abd = real array of dimension abd(nabd,n). ! on return contains the values of the matrix stored in ! banded form. The j-th column of abd contains the elements ! of the j-th column of the original matrix comprised in the ! band ( i in (j-ml,j+mu) ) with the lowest diagonal at ! the bottom row (row lowd). See details below for this format. ! ! ml = integer. equal to the bandwidth of the strict lower part of A ! mu = integer. equal to the bandwidth of the strict upper part of A ! if job=1 on entry then these two values are internally computed. ! ! lowd = integer. row number in abd where the lowest diagonal ! (leftmost) of A is located on return. In case lowd = 0 ! on return, then it is defined to ml+mu+1 on return and the ! lowd will contain this value on return. ` ! ! ierr = integer. used for error messages. On return: ! ierr == 0 :means normal return ! ierr == -1 : means invalid value for lowd. (either < 0 ! or larger than nabd). ! ierr == -2 : means that lowd is not large enough and as ! result the matrix cannot be stored in array abd. ! lowd should be at least ml+mu+1, where ml and mu are as ! provided on output. ! implicit none integer n integer nabd real ( kind = 8 ) a(*) real ( kind = 8 ) abd(nabd,n) integer i integer ia(n+1) integer ierr integer ii integer j integer ja(*) integer job integer k integer lowd integer m integer mdiag integer ml integer mu ! ! Determine ML and MU. ! ierr = 0 if ( job == 1 ) then call getbwd ( n, a, ja, ia, ml, mu ) end if m = ml + mu + 1 if ( lowd == 0 ) then lowd = m end if if ( lowd < m ) then ierr = -2 end if if ( nabd < lowd .or. lowd < 0 ) then ierr = -1 end if if ( ierr < 0 ) then return end if do i = 1, m ii = lowd - i + 1 abd(ii,1:n) = 0.0D+00 end do mdiag = lowd - ml do i = 1, n do k = ia(i), ia(i+1)-1 j = ja(k) abd(i-j+mdiag,j) = a(k) end do end do return end subroutine csrbsr ( n, nblk, na, a, ja, ia, ao, jao, iao ) !*****************************************************************************80 ! !! CSRBSR converts Compressed Sparse Row to Block Sparse Row. ! ! Discussion: ! ! This routine does the reverse of BSRCSR. It converts ! a matrix stored in a general compressed a, ja, ia format into a ! a block reduced matrix a(*,*),ja(*),ia(*) format. The code ! assumes that the original matrix is indeed a block format ! and that the elements are ordered in such a way that their ! column numbers are increasing. (This can be achieved ! by transposing a, ja, ia twice, putting the resulting matrix ! into a, ja, ia). ! ! See routine bsrcsr for more details on data structure for blocked ! matrices. The input matrix is a, ja, ia (in compressed format) and ! the output matrix is the matrix ao, jao, iao in block-reduced ! format. ! ! This code is not in place. ! ! See routine bsrcsr for details on data sctructure for ! block sparse row format. ! ! The routine assumes that the input matrix has been ! sorted in such a way that the column indices are always ! in increasing order for the same row. ! for all k "in the SAME ROW." ! ! THERE IS NO CHECKING AS TO WHETHER the input is correct. ! it is recommended to use the routine blchk to check ! if the matrix is a block-matrix before calling csrbsr. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the row dimension of the matrix. ! ! Input, integer NBLK, the dimension of each block. ! NBLK must divide N. ! ! na = first dimension of array ao as declared in calling program ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! on return: ! ! ao = real array containing the values of the matrix. For details ! on the format see below. Each row of a contains the nblk x nblk ! block matrix unpacked column-wise (this allows the user to ! declare the array a as a(na,nblk,nblk) on entry if desired). ! the block rows are stored in sequence just as for the compressed ! sparse row format. ! jao = integer array of length n/nblk. ja(k) contains the column index ! of the leading element, i.e., the element (1,1) of the block ! that is held in the row a(k,*) of the value array. ! iao = integer array of length n/nblk+1. ia(i) points to the beginning ! of block row number i in the arrays a and ja. ! implicit none integer n integer na real ( kind = 8 ) a(*) real ( kind = 8 ) ao(na,*) integer i integer i1 integer ia(n+1) integer iao(*) integer ibrow integer ii integer irow integer j integer j1 integer j2 integer ja(*) integer jao(*) integer k integer k1 integer len integer lena integer nblk integer nr ! ! NR is the dimension of the reduced matrix. ! nr = n / nblk iao(1) = 1 ibrow = 1 irow = 1 ! ! The main loop. ! do ii = 1, nr ! ! I1 = starting position for group of nblk rows in original matrix. ! i1 = ia(irow) ! ! LENA = length of each row in that group in the original matrix. ! lena = ia(irow+1) - i1 ! ! LEN = length of each block-row in that group in the output matrix. ! len = lena / nblk k1 = iao(ibrow) ! ! Copy the real values of A. ! ! For each block. ! do k = 0, len-1 ! ! Store column positions of the (1,1) elements of each block. ! jao(k1+k) = ja(i1+nblk*k) ! ! For each column. ! do j = 1, nblk j1 = ( j - 1 ) * nblk j2 = i1 + k * nblk + j - 1 ! ! For each row. ! do i = 1, nblk ao(k1+k,j1+i) = a(j2+(i-1)*lena) end do end do end do ! ! Done with a whole block row. ! Update IAO, IBROW and IROW. ! iao(ibrow+1) = iao(ibrow) + len ibrow = ibrow + 1 irow = irow + nblk end do return end subroutine csrcoo ( nrow, job, nzmax, a, ja, ia, nnz, ao, ir, jc, ierr ) !*****************************************************************************80 ! !! CSRCOO converts Compressed Sparse Row to Coordinate format. ! ! Discussion: ! ! This routine converts a matrix that is stored in row general sparse ! A, JA, IA format into coordinate format AO, IR, JC. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! job = integer serving as a job indicator. ! if job = 1 fill in only the array ir, ignore jc, and ao. ! if job = 2 fill in ir, and jc but not ao ! if job = 3 fill in everything. ! The reason why these options are provided is that on return ! ao and jc are the same as a, ja. So when job = 3, a and ja are ! simply copied into ao, jc. When job=2, only jc and ir are ! returned. With job=1 only the array ir is returned. Moreover, ! the algorithm is in place: ! call csrcoo (nrow,1,nzmax,a,ja,ia,nnz,a,ia,ja,ierr) ! will write the output matrix in coordinate format on a, ja,ia. ! (Important: note the order in the output arrays a, ja, ia. ) ! i.e., ao can be the same as a, ir can be the same as ia ! and jc can be the same as ja. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! nzmax = length of space available in ao, ir, jc. ! the code will stop immediatly if the number of ! nonzero elements found in input matrix exceeds nzmax. ! ! on return: !- ! ao, ir, jc = matrix in coordinate format. ! ! nnz = number of nonzero elements in matrix. ! ! ierr = integer error indicator. ! ierr == 0 means normal retur ! ierr == 1 means that the the code stopped ! because there was no space in ao, ir, jc ! (according to the value of nzmax). ! implicit none integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) ao(*) integer i integer ia(nrow+1) integer ierr integer ir(*) integer ja(*) integer jc(*) integer job integer k integer k1 integer k2 integer nnz integer nzmax ierr = 0 nnz = ia(nrow+1)-1 if ( nzmax < nnz ) then ierr = 1 return end if if ( 3 <= job ) then ao(1:nnz) = a(1:nnz) end if if ( 2 <= job ) then jc(1:nnz) = ja(1:nnz) end if ! ! Copy backward. ! do i = nrow, 1, -1 k1 = ia(i+1) - 1 k2 = ia(i) do k = k1, k2, -1 ir(k) = i end do end do return end subroutine csrcsc ( n, job, ipos, a, ja, ia, ao, jao, iao ) !*****************************************************************************80 ! !! CSRCSC converts Compressed Sparse Row to Compressed Sparse Column. ! ! Discussion: ! ! This is essentially a transposition operation. ! ! It is NOT an in-place algorithm. ! ! This routine transposes a matrix stored in a, ja, ia format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input, integer JOB, indicates whether or not to fill the values of the ! matrix AO or only the pattern (IA, and JA). Enter 1 for yes. ! ! ipos = starting position in ao, jao of the transposed matrix. ! the iao array takes this into account (thus iao(1) is set to ipos.) ! Note: this may be useful if one needs to append the data structure ! of the transpose to that of A. In this case use ! call csrcsc (n,1,n+2,a,ja,ia,a,ja,ia(n+2)) ! for any other normal usage, enter ipos=1. ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real AO(*), JAO(*), IAO(N+1), the matrix in CSC ! Compressed Sparse Column format. ! implicit none integer n real ( kind = 8 ) a(*) real ( kind = 8 ) ao(*) integer i integer ia(n+1) integer iao(n+1) integer ipos integer j integer ja(*) integer jao(*) integer job integer k integer next ! ! Compute lengths of rows of A'. ! iao(1:n+1) = 0 do i = 1, n do k = ia(i), ia(i+1)-1 j = ja(k) + 1 iao(j) = iao(j) + 1 end do end do ! ! Compute pointers from lengths. ! iao(1) = ipos do i = 1, n iao(i+1) = iao(i) + iao(i+1) end do ! ! Do the actual copying. ! do i = 1, n do k = ia(i), ia(i+1)-1 j = ja(k) next = iao(j) if ( job == 1 ) then ao(next) = a(k) end if jao(next) = i iao(j) = next + 1 end do end do ! ! Reshift IAO and leave. ! do i = n, 1, -1 iao(i+1) = iao(i) end do iao(1) = ipos return end subroutine csrdia ( n, idiag, job, a, ja, ia, ndiag, diag, ioff, ao, & jao, iao, ind ) !*****************************************************************************80 ! !! CSRDIA converts Compressed Sparse Row to diagonal format. ! ! Discussion: ! ! This routine extracts IDIAG diagonals from the input matrix A, ! JA, IA, and puts the rest of the matrix in the output matrix AO, ! JAO, IAO. The diagonals to be extracted depend on the value of JOB. ! ! In the first case, the diagonals to be ! extracted are simply identified by their offsets provided in ioff ! by the caller. In the second case, the code internally determines ! the idiag most significant diagonals, i.e., those diagonals of the ! matrix which have the largest number of nonzero elements, and ! extracts them. ! ! The algorithm is in place: ao, jao, iao can be overwritten on ! a, ja, ia if desired. ! ! When the code is required to select the diagonals (job >= 10) ! the selection of the diagonals is done from left to right ! as a result if several diagonals have the same weight (number ! of nonzero elemnts) the leftmost one is selected first. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input/output, integer IDIAG. On intput, the number of diagonals ! to be extracted. On output, IDIAG may be modified to the ! actual number of diagonals found. ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! job = integer. serves as a job indicator. Job is better thought ! of as a two-digit number job=xy. If the first (x) digit ! is one on entry then the diagonals to be extracted are ! internally determined. In this case csrdia exctracts the ! idiag most important diagonals, i.e. those having the largest ! number on nonzero elements. If the first digit is zero ! then csrdia assumes that ioff(*) contains the offsets ! of the diagonals to be extracted. there is no verification ! that ioff(*) contains valid entries. ! The second (y) digit of job determines whether or not ! the remainder of the matrix is to be written on ao,jao,iao. ! If it is zero then ao, jao, iao is not filled, i.e., ! the diagonals are found and put in array diag and the rest is ! is discarded. if it is one, ao, jao, iao contains matrix ! of the remaining elements. ! Thus: ! job= 0 means do not select diagonals internally (pick those ! defined by ioff) and do not fill ao,jao,iao ! job= 1 means do not select diagonals internally ! and fill ao,jao,iao ! job=10 means select diagonals internally ! and do not fill ao,jao,iao ! job=11 means select diagonals internally ! and fill ao,jao,iao ! ! Input, integer NDIAG, the first dimension of array DIAG. ! ! on return: ! ! diag = real array of size (ndiag x idiag) containing the diagonals ! of A on return ! ! ioff = integer array of length idiag, containing the offsets of the ! diagonals to be extracted. ! ao, jao ! iao = remainder of the matrix in a, ja, ia format. ! ! work arrays: ! ! ind = integer array of length 2*n-1 used as integer work space. ! needed only when job>=10 i.e., in case the diagonals are to ! be selected internally. ! implicit none integer idiag integer ndiag real ( kind = 8 ) a(*) real ( kind = 8 ) ao(*) real ( kind = 8 ) diag(ndiag,idiag) integer i integer ia(*) integer iao(*) integer idum integer ii integer ind(*) integer ioff(*) integer j integer ja(*) integer jao(*) integer jmax integer job integer job1 integer job2 integer k integer ko integer l integer n integer n2 job1 = job / 10 job2 = job - job1 * 10 if ( job1 /= 0 ) then n2 = n + n - 1 call infdia ( n, ja, ia, ind, idum ) ! ! Determine the diagonals to accept. ! ii = 0 do ii = ii + 1 jmax = 0 do k = 1, n2 j = ind(k) if ( jmax < j ) then i = k jmax = j end if end do if ( jmax <= 0 ) then ii = ii - 1 exit end if ioff(ii) = i - n ind(i) = - jmax if ( idiag <= ii ) then exit end if end do idiag = ii end if ! ! Initialize DIAG to zero. ! diag(1:n,1:idiag) = 0.0D+00 ko = 1 ! ! Extract diagonals and accumulate remaining matrix. ! do i = 1, n do k = ia(i), ia(i+1)-1 j = ja(k) do l = 1, idiag if ( j - i == ioff(l) ) then diag(i,l) = a(k) go to 51 end if end do ! ! Append element not in any diagonal to AO, JAO, IAO. ! if ( job2 /= 0 ) then ao(ko) = a(k) jao(ko) = j ko = ko + 1 end if 51 continue end do if ( job2 /= 0 ) then ind(i+1) = ko end if end do ! ! Finish with IAO. ! if ( job2 /= 0 ) then iao(1) = 1 iao(2:n+1) = ind(2:n+1) end if return end subroutine csrdns ( nrow, ncol, a, ja, ia, dns, ndns, ierr ) !*****************************************************************************80 ! !! CSRDNS converts Compressed Sparse Row to Dense format. ! ! Discussion: ! ! This routine converts a row-stored sparse matrix into a densely stored one. ! ! Modified: ! ! 07 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, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Output, real DNS(NDNS,NDNS), the dense array containing a ! copy of the matrix. ! ! Input, integer NDNS, the dimension of the DNS array. ! ! Output, integer IERR, error indicator. ! 0, means normal return ! i, means that the code has stopped when processing ! row number i, because it found a column number > ncol. ! implicit none integer ncol integer ndns real ( kind = 8 ) a(*) real ( kind = 8 ) dns(ndns,ncol) integer i integer ia(*) integer ierr integer j integer ja(*) integer k integer nrow ierr = 0 dns(1:nrow,1:ncol) = 0.0D+00 do i = 1, nrow do k = ia(i), ia(i+1)-1 j = ja(k) if ( ncol < j ) then ierr = i return end if dns(i,j) = a(k) end do end do return end subroutine csrell ( nrow, a, ja, ia, maxcol, coef, jcoef, ncoef, & ndiag, ierr ) !*****************************************************************************80 ! !! CSRELL converts Compressed Sparse Row to Ellpack/Itpack format ! ! Discussion: ! ! This routine converts a matrix stored in the general A, JA, IA ! format into the COEF, JCOEF Ellpack/Itpack format. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix A. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix, stored in ! compressed sparse row format. ! ! Input, integer NCOEF, the first dimension of arrays COEF, and JCOEF. ! ! Input, integer MAXCOL, the number of columns available in COEF. ! ! Output, real COEF(NCOEF,MAXCOL), the values of the matrix A in ! Ellpack/Itpack format. ! ! Output, integer JCOEF(NCOEF,MAXCOL), the column indices of each entry ! in COEF. ! ! Output, integer NDIAG, the number of active 'diagonals' found. ! ! Output, integer IERR, an error flag. ! 0 = correct return. ! nonzero means that NDIAG, the number of diagonals found, exceeds ! the limit of MAXCOL. ! implicit none integer maxcol integer ncoef integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) coef(ncoef,maxcol) integer i integer ia(nrow+1) integer ierr integer ja(*) integer jcoef(ncoef,maxcol) integer k integer k1 integer k2 integer ndiag ! ! Determine the length of each row of the lower part of A. ! ierr = 0 ndiag = 0 do i = 1, nrow k = ia(i+1) - ia(i) ndiag = max ( ndiag, k ) end do ! ! Check that sufficient columns are available. ! if ( maxcol < ndiag ) then ierr = 1 return end if ! ! Initialize COEF and JCOEF. ! coef(1:nrow,1:ndiag) = 0.0D+00 jcoef(1:nrow,1:ndiag) = i ! ! Copy elements by row. ! do i = 1, nrow k1 = ia(i) k2 = ia(i+1)-1 do k = k1, k2 coef(i,k-k1+1) = a(k) jcoef(i,k-k1+1) = ja(k) end do end do return end subroutine csrjad ( nrow, a, ja, ia, idiag, iperm, ao, jao, iao ) !*****************************************************************************80 ! !! CSRJAD converts Compressed Sparse Row to Jagged Diagonal storage. ! ! Discussion: ! ! This routine converts a matrix stored in the compressed sparse ! row format to the jagged diagonal format. The data structure ! for the JAD (Jagged Diagonal storage) is as follows. The rows of ! the matrix are implicitly permuted so that their lengths are in ! decreasing order. The real entries AO(*) and their column indices ! JAO(*) are stored in succession. The number of such diagonals is IDIAG. ! The lengths of each of these diagonals is stored in IAO(*). ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Reference: ! ! E. Anderson, Youcef Saad, ! Solving sparse triangular systems on parallel computers, ! International Journal of High Speed Computing, ! Volume 1, pages 73-96, 1989. ! ! Youcef Saad, ! Krylov Subspace Methods on Supercomputers, ! SIAM Journal on Statistical and Scientific Computing, ! Volume 10, pages 1200-1232, 1989. ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix in CSR ! Compressed Sparse Row format. ! ! on return: ! ! Output, integer IDIAG, the number of jagged diagonals in the data ! structure A, JA, IA. ! ! Output, integer IPERM(NROW), the permutation of the rows that leads ! to a decreasing order of the number of nonzero elements. ! ! ao = real array containing the values of the matrix A in ! jagged diagonal storage. The j-diagonals are stored ! in ao in sequence. ! ! Output, integer JAO(*), the column indices of the entries in ao. ! ! iao = integer array containing pointers to the beginning ! of each j-diagonal in ao, jao. iao is also used as ! a work array and it should be of length n at least. ! implicit none integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) ao(*) integer i integer ia(nrow+1) integer iao(nrow) integer idiag integer ilo integer iperm(nrow) integer j integer ja(*) integer jao(*) integer jj integer k integer k0 integer k1 integer len ! ! Define initial IPERM and get lengths of each row. ! JAO is used a work vector to store tehse lengths. ! idiag = 0 ilo = nrow do j = 1, nrow iperm(j) = j len = ia(j+1) - ia(j) ilo = min ( ilo, len ) idiag = max ( idiag, len ) jao(j) = len end do ! ! Call the sorter to get permutation. Use IAO as a work array. ! call dcsort ( jao, nrow, iao, iperm, ilo, idiag ) ! ! Define the output data structure. First lengths of the J-diagonals. ! iao(1:nrow) = 0 do k = 1, nrow len = jao(iperm(k)) do i = 1, len iao(i) = iao(i) + 1 end do end do ! ! Get the output matrix itself. ! k1 = 1 k0 = k1 do jj = 1, idiag len = iao(jj) do k = 1, len i = ia(iperm(k)) + jj -1 ao(k1) = a(i) jao(k1) = ja(i) k1 = k1 + 1 end do iao(jj) = k0 k0 = k1 end do iao(idiag+1) = k1 return end subroutine csrlnk ( n, a, ja, ia, link ) !*****************************************************************************80 ! !! CSRLNK converts Compressed Sparse Row to Linked storage format. ! ! Discussion: ! ! This routine translates a matrix stored in compressed sparse ! row into one with a linked list storage format. Only the link ! array needs to be obtained since the arrays A, JA, and IA may ! be unchanged and have carry the same meaning for the output matrix. ! ! In other words a, ja, ia, link ia the output linked list data ! structure with a, ja, ia being the same. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! on return: ! ! a = nonzero elements. ! ! ja = column positions. ! ! ia = points to the first row of matrix in structure. ! ! link = integer array of size containing the linked list information. ! link(k) points to the next element of the row after element ! ao(k), jcol(k). if link(k) = 0, then there is no next element, ! i.e., ao(k), jcol(k) is the last element of the current row. ! implicit none integer n real ( kind = 8 ) a(*) integer i integer ia(n+1) integer ja(*) integer k integer link(*) ! ! Loop through all rows. ! do i = 1, n do k = ia(i), ia(i+1)-2 link(k) = k + 1 end do link(ia(i+1)-1) = 0 end do return end subroutine csrmsr ( n, a, ja, ia, ao, jao, wk, iwk ) !*****************************************************************************80 ! !! CSRMSR converts Compressed Sparse Row to Modified Sparse Row. ! ! Discussion: ! ! This routine converts a general sparse matrix a, ja, ia into ! a compressed matrix using a separated diagonal (referred to as ! the bell-labs format as it is used by bell labs semi conductor ! group. We refer to it here as the modified sparse row format. ! ! This has been coded in such a way that one can overwrite ! the output matrix onto the input matrix if desired by a call of ! the form ! ! call csrmsr (n, a, ja, ia, a, ja, wk,iwk) ! ! In case ao, jao, are different from a, ja, then one can ! use ao, jao as the work arrays in the calling sequence: ! ! call csrmsr (n, a, ja, ia, ao, jao, ao,jao) ! ! Algorithm is in place. i.e. both: ! ! call csrmsr (n, a, ja, ia, ao, jao, ao,jao) ! (in which ao, jao, are different from a, ja) ! and ! call csrmsr (n, a, ja, ia, a, ja, wk,iwk) ! (in which wk, jwk, are different from a, ja) ! are OK. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! on return : ! ! ao, jao = sparse matrix in modified sparse row storage format: ! + ao(1:n) contains the diagonal of the matrix. ! + ao(n+2:nnz) contains the nondiagonal elements of the ! matrix, stored rowwise. ! + jao(n+2:nnz) : their column indices ! + jao(1:n+1) contains the pointer array for the nondiagonal ! elements in ao(n+1:nnz) and jao(n+2:nnz). ! i.e., for i <= n+1 jao(i) points to beginning of row i ! in arrays ao, jao. ! here nnz = number of nonzero elements+1 ! ! Work array, real WK(N). ! ! Work array, integer IWK(N+1). ! implicit none integer n real ( kind = 8 ) a(*) real ( kind = 8 ) ao(*) integer i integer ia(n+1) integer icount integer ii integer iptr integer iwk(n+1) integer j integer ja(*) integer jao(*) integer k real ( kind = 8 ) wk(n) icount = 0 ! ! Store away diagonal elements and count nonzero diagonal elements. ! do i = 1, n wk(i) = 0.0D+00 iwk(i+1) = ia(i+1) - ia(i) do k = ia(i), ia(i+1)-1 if ( ja(k) == i ) then wk(i) = a(k) icount = icount + 1 iwk(i+1) = iwk(i+1) - 1 end if end do end do ! ! Compute total length. ! iptr = n + ia(n+1) - icount ! ! Copy backwards, to avoid collisions. ! do ii = n, 1, -1 do k = ia(ii+1)-1, ia(ii), -1 j = ja(k) if ( j /= ii ) then ao(iptr) = a(k) jao(iptr) = j iptr = iptr - 1 end if end do end do ! ! Compute the pointer values and copy WK. ! jao(1) = n + 2 do i = 1, n ao(i) = wk(i) jao(i+1) = jao(i) + iwk(i+1) end do return end subroutine csrncf ( nrow, a, ja, ia, maxnz, nonz, coef, jcoef, ierr ) !*****************************************************************************80 ! !! CSRNCF converts CSR to NSPCG NCF format. ! ! Discussion: ! ! This routine converts a matrix stored in the general A, JA, IA ! compressed sparse row format into the Nonsymmetric Coordinate Format ! used as storage format 5 by NSPCG. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NROW, the row dimension of the matrix A. ! ! Input, real A(*), integer JA(*), IA(NROW+1), the matrix, stored in ! compressed sparse row format. ! ! Input, integer MAXNZ, the maximum number of nonzeros allowed for ! in the storage of COEF and JCOEF. ! ! Output, integer NONZ, the actual number of nonzeros encountered. ! ! Output, real COEF(MAXNZ), the values of the matrix A in NCF format. ! ! Output, integer JCOEF(MAXNZ,2), the row and column indices of each ! entry in COEF. ! ! Output, integer IERR, an error flag. ! 0 = correct return. ! nonzero means that MAXNZ < NONZ. ! implicit none integer maxnz integer nrow real ( kind = 8 ) a(*) real ( kind = 8 ) coef(maxnz) integer i integer ia(nrow+1) integer ierr integer ja(*) integer jcoef(maxnz,2) integer k integer k1 integer k2 integer nonz ierr = 0 ! ! Initialize COEF and JCOEF. ! coef(1:maxnz) = 0.0D+00 jcoef(1:maxnz,1:2) = 0 ! ! The first N entries are reserved for the diagonals. ! do i = 1, nrow jcoef(i,1:2) = i end do nonz = nrow if ( maxnz < nonz ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CSRNCF - Fatal error!' write ( *, '(a)' ) ' MAXNZ < NONZ.' ierr = 1 return end if do i = 1, nrow k1 = ia(i) k2 = ia(i+1) - 1 do k = k1, k2 if ( ja(k) == i ) then coef(i) = coef(i) + a(k) else if ( 0.0D+00 /= a(k) ) then nonz = nonz + 1 if ( maxnz < nonz ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CSRNCF - Fatal error!' write ( *, '(a)' ) ' MAXNZ < NONZ.' ierr = 1 return end if coef(nonz) = a(k) jcoef(nonz,1) = i jcoef(nonz,2) = ja(k) end if end do end do return end subroutine csrssk ( n, imod, a, ja, ia, asky, isky, nzmax, ierr ) !*****************************************************************************80 ! !! CSRSSK converts Compressed Sparse Row to Symmetric Skyline Format. ! ! Discussion: ! ! This routine translates a compressed sparse row or a symmetric ! sparse row format into a symmetric skyline format. ! the input matrix can be in either compressed sparse row or the ! symmetric sparse row format. The output matrix is in a symmetric ! skyline format: a real array containing the (active portions) of the ! rows in sequence and a pointer to the beginning of each row. ! ! This module is NOT in place. ! ! Even when imod = 2, length of isky is n+1, not n. ! ! Modified: ! ! 07 January 2004 ! ! Author: ! ! Youcef Saad ! ! Parameters: ! ! Input, integer N, the order of the matrix. ! ! imod = integer indicating the variant of skyline format wanted: ! imod = 0 means the pointer isky points to the `zeroth' ! element of the row, i.e., to the position of the diagonal ! element of previous row (for i = 1, isky(1)= 0) ! imod = 1 means that itpr points to the beginning of the row. ! imod = 2 means that isky points to the end of the row (diagonal ! element) ! ! Input, real A(*), integer JA(*), IA(N+1), the matrix in CSR ! Compressed Sparse Row format. ! ! Input, integer NZMAX, the amount of storage available in ASKY. ! ! on return: ! ! asky = real array containing the values of the matrix stored in skyline ! format. asky contains the sequence of active rows from ! i = 1, to n, an active row being the row of elemnts of ! the matrix contained between the leftmost nonzero element ! and the diagonal element. ! ! isky = integer array of size n+1 containing the pointer array to ! each row. The meaning of isky depends on the input value of ! imod (see above). ! ! ierr = integer. Error message