c----------------------------------------------------------------------c c S P A R S K I T c c----------------------------------------------------------------------c c FORMAT CONVERSION MODULE c c----------------------------------------------------------------------c c contents: c c---------- c c csrdns : converts a row-stored sparse matrix into the dense format. c c dnscsr : converts a dense matrix to a sparse storage format. c c coocsr : converts coordinate to to csr format c c coicsr : in-place conversion of coordinate to csr format c c csrcoo : converts compressed sparse row to coordinate. c c csrssr : converts compressed sparse row to symmetric sparse row c c ssrcsr : converts symmetric sparse row to compressed sparse row c c csrell : converts compressed sparse row to ellpack format c c ellcsr : converts ellpack format to compressed sparse row format c c csrmsr : converts compressed sparse row format to modified sparse c c row format c c msrcsr : converts modified sparse row format to compressed sparse c c row format. c c csrcsc : converts compressed sparse row format to compressed sparse c c column format (transposition) c c csrcsc2 : rectangular version of csrcsc c c csrlnk : converts compressed sparse row to linked list format c c lnkcsr : converts linked list format to compressed sparse row fmt c c csrdia : converts a compressed sparse row format into a diagonal c c format. c c diacsr : converts a diagonal format into a compressed sparse row c c format. c c bsrcsr : converts a block-row sparse format into a compressed c c sparse row format. c c csrbsr : converts a compressed sparse row format into a block-row c c sparse format. c c csrbnd : converts a compressed sparse row format into a banded c c format (linpack style). c c bndcsr : converts a banded format (linpack style) into a compressed c c sparse row storage. c c csrssk : converts the compressed sparse row format to the symmetric c c skyline format c c sskssr : converts symmetric skyline format to symmetric sparse row c c format. c c csrjad : converts the csr format into the jagged diagonal format c c jadcsr : converts the jagged-diagonal format into the csr format c c csruss : Compressed Sparse Row to Unsymmetric Sparse Skyline c c format c c usscsr : Unsymmetric Sparse Skyline format to Compressed Sparse Row c c csrsss : Compressed Sparse Row to Symmetric Sparse Skyline format c c ssscsr : Symmetric Sparse Skyline format to Compressed Sparse Row c c csrvbr : Converts compressed sparse row to var block row format c c vbrcsr : Converts var block row to compressed sparse row format c c csorted : Checks if matrix in CSR format is sorted by columns c c--------- miscalleneous additions not involving the csr format--------c c cooell : converts coordinate to Ellpack/Itpack format c c dcsort : sorting routine used by crsjad c c----------------------------------------------------------------------c subroutine csrdns(nrow,ncol,a,ja,ia,dns,ndns,ierr) real*8 dns(ndns,*),a(*) integer ja(*),ia(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Dense c----------------------------------------------------------------------- c c converts a row-stored sparse matrix into a densely stored one c c On entry: c---------- c c nrow = row-dimension of a c ncol = column dimension of a c a, c ja, c ia = input matrix in compressed sparse row format. c (a=value array, ja=column array, ia=pointer array) c dns = array where to store dense matrix c ndns = first dimension of array dns c c on return: c----------- c dns = the sparse matrix a, ja, ia has been stored in dns(ndns,*) c c ierr = integer error indicator. c ierr .eq. 0 means normal return c ierr .eq. i means that the code has stopped when processing c row number i, because it found a column number .gt. ncol. c c----------------------------------------------------------------------- ierr = 0 do 1 i=1, nrow do 2 j=1,ncol dns(i,j) = 0.0d0 2 continue 1 continue c do 4 i=1,nrow do 3 k=ia(i),ia(i+1)-1 j = ja(k) if (j .gt. ncol) then ierr = i return endif dns(i,j) = a(k) 3 continue 4 continue return c---- end of csrdns ---------------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine dnscsr(nrow,ncol,nzmax,dns,ndns,a,ja,ia,ierr) real*8 dns(ndns,*),a(*) integer ia(*),ja(*) c----------------------------------------------------------------------- c Dense to Compressed Row Sparse c----------------------------------------------------------------------- c c converts a densely stored matrix into a row orientied c compactly sparse matrix. ( reverse of csrdns ) c Note: this routine does not check whether an element c is small. It considers that a(i,j) is zero if it is exactly c equal to zero: see test below. c----------------------------------------------------------------------- c on entry: c--------- c c nrow = row-dimension of a c ncol = column dimension of a c nzmax = maximum number of nonzero elements allowed. This c should be set to be the lengths of the arrays a and ja. c dns = input nrow x ncol (dense) matrix. c ndns = first dimension of dns. c c on return: c---------- c c a, ja, ia = value, column, pointer arrays for output matrix c c ierr = integer error indicator: c ierr .eq. 0 means normal retur c ierr .eq. i means that the the code stopped while c processing row number i, because there was no space left in c a, and ja (as defined by parameter nzmax). c----------------------------------------------------------------------- ierr = 0 next = 1 ia(1) = 1 do 4 i=1,nrow do 3 j=1, ncol if (dns(i,j) .eq. 0.0d0) goto 3 if (next .gt. nzmax) then ierr = i return endif ja(next) = j a(next) = dns(i,j) next = next+1 3 continue ia(i+1) = next 4 continue return c---- end of dnscsr ---------------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine coocsr(nrow,nnz,a,ir,jc,ao,jao,iao) c----------------------------------------------------------------------- real*8 a(*),ao(*),x integer ir(*),jc(*),jao(*),iao(*) c----------------------------------------------------------------------- c Coordinate to Compressed Sparse Row c----------------------------------------------------------------------- c converts a matrix that is stored in coordinate format c a, ir, jc into a row general sparse ao, jao, iao format. c c on entry: c--------- c nrow = dimension of the matrix c nnz = number of nonzero elements in matrix c a, c ir, c jc = matrix in coordinate format. a(k), ir(k), jc(k) store the nnz c nonzero elements of the matrix with a(k) = actual real value of c the elements, ir(k) = its row number and jc(k) = its column c number. The order of the elements is arbitrary. c c on return: c----------- c ir is destroyed c c ao, jao, iao = matrix in general sparse matrix format with ao c continung the real values, jao containing the column indices, c and iao being the pointer to the beginning of the row, c in arrays ao, jao. c c Notes: c------ This routine is NOT in place. See coicsr c c------------------------------------------------------------------------ do 1 k=1,nrow+1 iao(k) = 0 1 continue c determine row-lengths. do 2 k=1, nnz iao(ir(k)) = iao(ir(k))+1 2 continue c starting position of each row.. k = 1 do 3 j=1,nrow+1 k0 = iao(j) iao(j) = k k = k+k0 3 continue c go through the structure once more. Fill in output matrix. do 4 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 4 continue c shift back iao do 5 j=nrow,1,-1 iao(j+1) = iao(j) 5 continue iao(1) = 1 return c------------- end of coocsr ------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine coicsr (n,nnz,job,a,ja,ia,iwk) integer ia(nnz),ja(nnz),iwk(n+1) real*8 a(*) c------------------------------------------------------------------------ c IN-PLACE coo-csr conversion routine. c------------------------------------------------------------------------ c this subroutine converts a matrix stored in coordinate format into c the csr format. The conversion is done in place in that the arrays c a,ja,ia of the result are overwritten onto the original arrays. c------------------------------------------------------------------------ c on entry: c--------- c n = integer. row dimension of A. c nnz = integer. number of nonzero elements in A. c job = integer. Job indicator. when job=1, the real values in a are c filled. Otherwise a is not touched and the structure of the c array only (i.e. ja, ia) is obtained. c a = real array of size nnz (number of nonzero elements in A) c containing the nonzero elements c ja = integer array of length nnz containing the column positions c of the corresponding elements in a. c ia = integer array of length nnz containing the row positions c of the corresponding elements in a. c iwk = integer work array of length n+1 c on return: c---------- c a c ja c ia = contains the compressed sparse row data structure for the c resulting matrix. c Note: c------- c the entries of the output matrix are not sorted (the column c indices in each are not in increasing order) use coocsr c if you want them sorted. c----------------------------------------------------------------------c c Coded by Y. Saad, Sep. 26 1989 c c----------------------------------------------------------------------c real*8 t,tnext logical values c----------------------------------------------------------------------- values = (job .eq. 1) c find pointer array for resulting matrix. do 35 i=1,n+1 iwk(i) = 0 35 continue do 4 k=1,nnz i = ia(k) iwk(i+1) = iwk(i+1)+1 4 continue c------------------------------------------------------------------------ iwk(1) = 1 do 44 i=2,n iwk(i) = iwk(i-1) + iwk(i) 44 continue c c loop for a cycle in chasing process. c init = 1 k = 0 5 if (values) t = a(init) i = ia(init) j = ja(init) ia(init) = -1 c------------------------------------------------------------------------ 6 k = k+1 c current row number is i. determine where to go. ipos = iwk(i) c save the chased element. if (values) tnext = a(ipos) inext = ia(ipos) jnext = ja(ipos) c then occupy its location. if (values) a(ipos) = t ja(ipos) = j c update pointer information for next element to come in row i. iwk(i) = ipos+1 c determine next element to be chased, if (ia(ipos) .lt. 0) goto 65 t = tnext i = inext j = jnext ia(ipos) = -1 if (k .lt. nnz) goto 6 goto 70 65 init = init+1 if (init .gt. nnz) goto 70 if (ia(init) .lt. 0) goto 65 c restart chasing -- goto 5 70 do 80 i=1,n ia(i+1) = iwk(i) 80 continue ia(1) = 1 return c----------------- end of coicsr ---------------------------------------- c------------------------------------------------------------------------ end c----------------------------------------------------------------------- subroutine csrcoo (nrow,job,nzmax,a,ja,ia,nnz,ao,ir,jc,ierr) c----------------------------------------------------------------------- real*8 a(*),ao(*) integer ir(*),jc(*),ja(*),ia(nrow+1) c----------------------------------------------------------------------- c Compressed Sparse Row to Coordinate c----------------------------------------------------------------------- c converts a matrix that is stored in coordinate format c a, ir, jc into a row general sparse ao, jao, iao format. c c on entry: c--------- c nrow = dimension of the matrix. c job = integer serving as a job indicator. c if job = 1 fill in only the array ir, ignore jc, and ao. c if job = 2 fill in ir, and jc but not ao c if job = 3 fill in everything. c The reason why these options are provided is that on return c ao and jc are the same as a, ja. So when job = 3, a and ja are c simply copied into ao, jc. When job=2, only jc and ir are c returned. With job=1 only the array ir is returned. Moreover, c the algorithm is in place: c call csrcoo (nrow,1,nzmax,a,ja,ia,nnz,a,ia,ja,ierr) c will write the output matrix in coordinate format on a, ja,ia. c c a, c ja, c ia = matrix in compressed sparse row format. c nzmax = length of space available in ao, ir, jc. c the code will stop immediatly if the number of c nonzero elements found in input matrix exceeds nzmax. c c on return: c----------- c ao, ir, jc = matrix in coordinate format. c c nnz = number of nonzero elements in matrix. c ierr = integer error indicator. c ierr .eq. 0 means normal retur c ierr .eq. 1 means that the the code stopped c because there was no space in ao, ir, jc c (according to the value of nzmax). c c NOTES: 1)This routine is PARTIALLY in place: csrcoo can be called with c ao being the same array as as a, and jc the same array as ja. c but ir CANNOT be the same as ia. c 2) note the order in the output arrays, c------------------------------------------------------------------------ ierr = 0 nnz = ia(nrow+1)-1 if (nnz .gt. nzmax) then ierr = 1 return endif c------------------------------------------------------------------------ goto (3,2,1) job 1 do 10 k=1,nnz ao(k) = a(k) 10 continue 2 do 11 k=1,nnz jc(k) = ja(k) 11 continue c c copy backward to allow for in-place processing. c 3 do 13 i=nrow,1,-1 k1 = ia(i+1)-1 k2 = ia(i) do 12 k=k1,k2,-1 ir(k) = i 12 continue 13 continue return c------------- end-of-csrcoo ------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrssr (nrow,a,ja,ia,nzmax,ao,jao,iao,ierr) real*8 a(*), ao(*), t integer ia(*), ja(*), iao(*), jao(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Symmetric Sparse Row c----------------------------------------------------------------------- c this subroutine extracts the lower triangular part of a matrix. c It can used as a means for converting a symmetric matrix for c which all the entries are stored in sparse format into one c in which only the lower part is stored. The routine is in place in c that the output matrix ao, jao, iao can be overwritten on c the input matrix a, ja, ia if desired. Csrssr has been coded to c put the diagonal elements of the matrix in the last position in c each row (i.e. in position ao(ia(i+1)-1 of ao and jao) c----------------------------------------------------------------------- c On entry c----------- c nrow = dimension of the matrix a. c a, ja, c ia = matrix stored in compressed row sparse format c c nzmax = length of arrays ao, and jao. c c On return: c----------- c ao, jao, c iao = lower part of input matrix (a,ja,ia) stored in compressed sparse c row format format. c c ierr = integer error indicator. c ierr .eq. 0 means normal return c ierr .eq. i means that the code has stopped when processing c row number i, because there is not enough space in ao, jao c (according to the value of nzmax) c c----------------------------------------------------------------------- ierr = 0 ko = 0 c----------------------------------------------------------------------- do 7 i=1, nrow kold = ko kdiag = 0 do 71 k = ia(i), ia(i+1) -1 if (ja(k) .gt. i) goto 71 ko = ko+1 if (ko .gt. nzmax) then ierr = i return endif ao(ko) = a(k) jao(ko) = ja(k) if (ja(k) .eq. i) kdiag = ko 71 continue if (kdiag .eq. 0 .or. kdiag .eq. ko) goto 72 c c exchange c t = ao(kdiag) ao(kdiag) = ao(ko) ao(ko) = t c k = jao(kdiag) jao(kdiag) = jao(ko) jao(ko) = k 72 iao(i) = kold+1 7 continue c redefine iao(n+1) iao(nrow+1) = ko+1 return c--------- end of csrssr ----------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine ssrcsr(job, value2, nrow, a, ja, ia, nzmax, & ao, jao, iao, indu, iwk, ierr) c .. Scalar Arguments .. integer ierr, job, nrow, nzmax, value2 c .. c .. Array Arguments .. integer ia(nrow+1), iao(nrow+1), indu(nrow), & iwk(nrow+1), ja(*), jao(nzmax) real*8 a(*), ao(nzmax) c .. c----------------------------------------------------------------------- c Symmetric Sparse Row to Compressed Sparse Row format c----------------------------------------------------------------------- c This subroutine converts a given matrix in SSR format to regular c CSR format by computing Ao = A + A' - diag(A), where A' is A c transpose. c c Typically this routine is used to expand the SSR matrix of c Harwell Boeing matrices, or to obtain a symmetrized graph of c unsymmetric matrices. c c This routine is inplace, i.e., (Ao,jao,iao) may be same as c (a,ja,ia). c c It is possible to input an arbitrary CSR matrix to this routine, c since there is no syntactical difference between CSR and SSR c format. It also removes duplicate entries and perform a partial c ordering. The output matrix has an order of lower half, main c diagonal and upper half after the partial ordering. c----------------------------------------------------------------------- c on entry: c--------- c c job = options c 0 -- duplicate entries are not removed. If the input matrix is c SSR (not an arbitary CSR) matrix, no duplicate entry should c arise from this routine. c 1 -- eliminate duplicate entries, zero entries. c 2 -- eliminate duplicate entries and perform partial ordering. c 3 -- eliminate duplicate entries, sort the entries in the c increasing order of clumn indices. c c value2= will the values of A be copied? c 0 -- only expand the graph (a, ao are not touched) c 1 -- expand the matrix with the values. c c nrow = column dimension of inout matrix c a, c ia, c ja = matrix in compressed sparse row format. c c nzmax = size of arrays ao and jao. SSRCSR will abort if the storage c provided in ao, jao is not sufficient to store A. See ierr. c c on return: c---------- c ao, jao, iao c = output matrix in compressed sparse row format. The resulting c matrix is symmetric and is equal to A+A'-D. ao, jao, iao, c can be the same as a, ja, ia in the calling sequence. c c indu = integer array of length nrow. INDU will contain pointers c to the beginning of upper traigular part if job > 1. c Otherwise it is also used as a work array (size nrow). c c iwk = integer work space (size nrow+1). c c ierr = integer. Serving as error message. If the length of the arrays c ao, jao exceeds nzmax, ierr returns the minimum value c needed for nzmax. otherwise ierr=0 (normal return). c c----------------------------------------------------------------------- c .. Local Scalars .. integer i, ipos, j, k, kfirst, klast, ko, kosav, nnz real*8 tmp c .. c .. Executable Statements .. ierr = 0 do 10 i = 1, nrow indu(i) = 0 iwk(i) = 0 10 continue iwk(nrow+1) = 0 c c .. compute number of elements in each row of (A'-D) c put result in iwk(i+1) for row i. c do 30 i = 1, nrow do 20 k = ia(i), ia(i+1) - 1 j = ja(k) if (j.ne.i) & iwk(j+1) = iwk(j+1) + 1 20 continue 30 continue c c .. find addresses of first elements of ouput matrix. result in iwk c iwk(1) = 1 do 40 i = 1, nrow indu(i) = iwk(i) + ia(i+1) - ia(i) iwk(i+1) = iwk(i+1) + indu(i) indu(i) = indu(i) - 1 40 continue c.....Have we been given enough storage in ao, jao ? nnz = iwk(nrow+1) - 1 if (nnz.gt.nzmax) then ierr = nnz return endif c c .. copy the existing matrix (backwards). c kosav = iwk(nrow+1) do 60 i = nrow, 1, -1 klast = ia(i+1) - 1 kfirst = ia(i) iao(i+1) = kosav kosav = iwk(i) ko = iwk(i) - kfirst iwk(i) = ko + klast + 1 do 50 k = klast, kfirst, -1 if (value2.ne.0) & ao(k+ko) = a(k) jao(k+ko) = ja(k) 50 continue 60 continue iao(1) = 1 c c now copy (A'-D). Go through the structure of ao, jao, iao c that has already been copied. iwk(i) is the address c of the next free location in row i for ao, jao. c do 80 i = 1, nrow do 70 k = iao(i), indu(i) j = jao(k) if (j.ne.i) then ipos = iwk(j) if (value2.ne.0) & ao(ipos) = ao(k) jao(ipos) = i iwk(j) = ipos + 1 endif 70 continue 80 continue if (job.le.0) return c c .. eliminate duplicate entries -- c array INDU is used as marker for existing indices, it is also the c location of the entry. c IWK is used to stored the old IAO array. c matrix is copied to squeeze out the space taken by the duplicated c entries. c do 90 i = 1, nrow indu(i) = 0 iwk(i) = iao(i) 90 continue iwk(nrow+1) = iao(nrow+1) k = 1 do 120 i = 1, nrow iao(i) = k ipos = iwk(i) klast = iwk(i+1) 100 if (ipos.lt.klast) then j = jao(ipos) if (indu(j).eq.0) then c .. new entry .. if (value2.ne.0) then if (ao(ipos) .ne. 0.0D0) then indu(j) = k jao(k) = jao(ipos) ao(k) = ao(ipos) k = k + 1 endif else indu(j) = k jao(k) = jao(ipos) k = k + 1 endif else if (value2.ne.0) then c .. duplicate entry .. ao(indu(j)) = ao(indu(j)) + ao(ipos) endif ipos = ipos + 1 go to 100 endif c .. remove marks before working on the next row .. do 110 ipos = iao(i), k - 1 indu(jao(ipos)) = 0 110 continue 120 continue iao(nrow+1) = k if (job.le.1) return c c .. partial ordering .. c split the matrix into strict upper/lower triangular c parts, INDU points to the the beginning of the strict upper part. c do 140 i = 1, nrow klast = iao(i+1) - 1 kfirst = iao(i) 130 if (klast.gt.kfirst) then if (jao(klast).lt.i .and. jao(kfirst).ge.i) then c .. swap klast with kfirst .. j = jao(klast) jao(klast) = jao(kfirst) jao(kfirst) = j if (value2.ne.0) then tmp = ao(klast) ao(klast) = ao(kfirst) ao(kfirst) = tmp endif endif if (jao(klast).ge.i) & klast = klast - 1 if (jao(kfirst).lt.i) & kfirst = kfirst + 1 go to 130 endif c if (jao(klast).lt.i) then indu(i) = klast + 1 else indu(i) = klast endif 140 continue if (job.le.2) return c c .. order the entries according to column indices c bubble-sort is used c do 190 i = 1, nrow do 160 ipos = iao(i), indu(i)-1 do 150 j = indu(i)-1, ipos+1, -1 k = j - 1 if (jao(k).gt.jao(j)) then ko = jao(k) jao(k) = jao(j) jao(j) = ko if (value2.ne.0) then tmp = ao(k) ao(k) = ao(j) ao(j) = tmp endif endif 150 continue 160 continue do 180 ipos = indu(i), iao(i+1)-1 do 170 j = iao(i+1)-1, ipos+1, -1 k = j - 1 if (jao(k).gt.jao(j)) then ko = jao(k) jao(k) = jao(j) jao(j) = ko if (value2.ne.0) then tmp = ao(k) ao(k) = ao(j) ao(j) = tmp endif endif 170 continue 180 continue 190 continue c return c---- end of ssrcsr ---------------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine xssrcsr (nrow,a,ja,ia,nzmax,ao,jao,iao,indu,ierr) integer ia(nrow+1),iao(nrow+1),ja(*),jao(nzmax),indu(nrow+1) real*8 a(*),ao(nzmax) c----------------------------------------------------------------------- c Symmetric Sparse Row to (regular) Compressed Sparse Row c----------------------------------------------------------------------- c this subroutine converts a symmetric matrix in which only the lower c part is stored in compressed sparse row format, i.e., c a matrix stored in symmetric sparse format, into a fully stored matrix c i.e., a matrix where both the lower and upper parts are stored in c compressed sparse row format. the algorithm is in place (i.e. result c may be overwritten onto the input matrix a, ja, ia ----- ). c the output matrix delivered by ssrcsr is such that each row starts with c the elements of the lower part followed by those of the upper part. c----------------------------------------------------------------------- c on entry: c--------- c c nrow = row dimension of inout matrix c a, c ia, c ja = matrix in compressed sparse row format. This is assumed to be c a lower triangular matrix. c c nzmax = size of arrays ao and jao. ssrcsr will abort if the storage c provided in a, ja is not sufficient to store A. See ierr. c c on return: c---------- c ao, iao, c jao = output matrix in compressed sparse row format. The resulting c matrix is symmetric and is equal to A+A**T - D, if c A is the original lower triangular matrix. ao, jao, iao, c can be the same as a, ja, ia in the calling sequence. c c indu = integer array of length nrow+1. If the input matrix is such c that the last element in each row is its diagonal element then c on return, indu will contain the pointers to the diagonal c element in each row of the output matrix. Otherwise used as c work array. c ierr = integer. Serving as error message. If the length of the arrays c ao, jao exceeds nzmax, ierr returns the minimum value c needed for nzmax. otherwise ierr=0 (normal return). c c----------------------------------------------------------------------- ierr = 0 do 1 i=1,nrow+1 indu(i) = 0 1 continue c c compute number of elements in each row of strict upper part. c put result in indu(i+1) for row i. c do 3 i=1, nrow do 2 k=ia(i),ia(i+1)-1 j = ja(k) if (j .lt. i) indu(j+1) = indu(j+1)+1 2 continue 3 continue c----------- c find addresses of first elements of ouput matrix. result in indu c----------- indu(1) = 1 do 4 i=1,nrow lenrow = ia(i+1)-ia(i) indu(i+1) = indu(i) + indu(i+1) + lenrow 4 continue c--------------------- enough storage in a, ja ? -------- nnz = indu(nrow+1)-1 if (nnz .gt. nzmax) then ierr = nnz return endif c c now copy lower part (backwards). c kosav = indu(nrow+1) do 6 i=nrow,1,-1 klast = ia(i+1)-1 kfirst = ia(i) iao(i+1) = kosav ko = indu(i) kosav = ko do 5 k = kfirst, klast ao(ko) = a(k) jao(ko) = ja(k) ko = ko+1 5 continue indu(i) = ko 6 continue iao(1) = 1 c c now copy upper part. Go through the structure of ao, jao, iao c that has already been copied (lower part). indu(i) is the address c of the next free location in row i for ao, jao. c do 8 i=1,nrow c i-th row is now in ao, jao, iao structure -- lower half part do 9 k=iao(i), iao(i+1)-1 j = jao(k) if (j .ge. i) goto 8 ipos = indu(j) ao(ipos) = ao(k) jao(ipos) = i indu(j) = indu(j) + 1 9 continue 8 continue return c----- end of xssrcsr -------------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrell (nrow,a,ja,ia,maxcol,coef,jcoef,ncoef, * ndiag,ierr) integer ia(nrow+1), ja(*), jcoef(ncoef,1) real*8 a(*), coef(ncoef,1) c----------------------------------------------------------------------- c Compressed Sparse Row to Ellpack - Itpack format c----------------------------------------------------------------------- c this subroutine converts matrix stored in the general a, ja, ia c format into the coef, jcoef itpack format. c c----------------------------------------------------------------------- c on entry: c---------- c nrow = row dimension of the matrix A. c c a, c ia, c ja = input matrix in compressed sparse row format. c c ncoef = first dimension of arrays coef, and jcoef. c c maxcol = integer equal to the number of columns available in coef. c c on return: c---------- c coef = real array containing the values of the matrix A in c itpack-ellpack format. c jcoef = integer array containing the column indices of coef(i,j) c in A. c ndiag = number of active 'diagonals' found. c c ierr = error message. 0 = correct return. If ierr .ne. 0 on c return this means that the number of diagonals found c (ndiag) exceeds maxcol. c c----------------------------------------------------------------------- c first determine the length of each row of lower-part-of(A) ierr = 0 ndiag = 0 do 3 i=1, nrow k = ia(i+1)-ia(i) ndiag = max0(ndiag,k) 3 continue c----- check whether sufficient columns are available. ----------------- if (ndiag .gt. maxcol) then ierr = 1 return endif c c fill coef with zero elements and jcoef with row numbers.------------ c do 4 j=1,ndiag do 41 i=1,nrow coef(i,j) = 0.0d0 jcoef(i,j) = i 41 continue 4 continue c c------- copy elements row by row.-------------------------------------- c do 6 i=1, nrow k1 = ia(i) k2 = ia(i+1)-1 do 5 k=k1,k2 coef(i,k-k1+1) = a(k) jcoef(i,k-k1+1) = ja(k) 5 continue 6 continue return c--- end of csrell------------------------------------------------------ c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine ellcsr(nrow,coef,jcoef,ncoef,ndiag,a,ja,ia,nzmax,ierr) integer ia(nrow+1), ja(*), jcoef(ncoef,1) real*8 a(*), coef(ncoef,1) c----------------------------------------------------------------------- c Ellpack - Itpack format to Compressed Sparse Row c----------------------------------------------------------------------- c this subroutine converts a matrix stored in ellpack-itpack format c coef-jcoef into the compressed sparse row format. It actually checks c whether an entry in the input matrix is a nonzero element before c putting it in the output matrix. The test does not account for small c values but only for exact zeros. c----------------------------------------------------------------------- c on entry: c---------- c c nrow = row dimension of the matrix A. c coef = array containing the values of the matrix A in ellpack format. c jcoef = integer arraycontains the column indices of coef(i,j) in A. c ncoef = first dimension of arrays coef, and jcoef. c ndiag = number of active columns in coef, jcoef. c c ndiag = on entry the number of columns made available in coef. c c on return: c---------- c a, ia, c ja = matrix in a, ia, ja format where. c c nzmax = size of arrays a and ja. ellcsr will abort if the storage c provided in a, ja is not sufficient to store A. See ierr. c c ierr = integer. serves are output error message. c ierr = 0 means normal return. c ierr = 1 means that there is not enough space in c a and ja to store output matrix. c----------------------------------------------------------------------- c first determine the length of each row of lower-part-of(A) ierr = 0 c-----check whether sufficient columns are available. ----------------- c c------- copy elements row by row.-------------------------------------- kpos = 1 ia(1) = kpos do 6 i=1, nrow do 5 k=1,ndiag if (coef(i,k) .ne. 0.0d0) then if (kpos .gt. nzmax) then ierr = kpos return endif a(kpos) = coef(i,k) ja(kpos) = jcoef(i,k) kpos = kpos+1 endif 5 continue ia(i+1) = kpos 6 continue return c--- end of ellcsr ----------------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrmsr (n,a,ja,ia,ao,jao,wk,iwk) real*8 a(*),ao(*),wk(n) integer ia(n+1),ja(*),jao(*),iwk(n+1) c----------------------------------------------------------------------- c Compressed Sparse Row to Modified - Sparse Row c Sparse row with separate main diagonal c----------------------------------------------------------------------- c converts a general sparse matrix a, ja, ia into c a compressed matrix using a separated diagonal (referred to as c the bell-labs format as it is used by bell labs semi conductor c group. We refer to it here as the modified sparse row format. c Note: this has been coded in such a way that one can overwrite c the output matrix onto the input matrix if desired by a call of c the form c c call csrmsr (n, a, ja, ia, a, ja, wk,iwk) c c In case ao, jao, are different from a, ja, then one can c use ao, jao as the work arrays in the calling sequence: c c call csrmsr (n, a, ja, ia, ao, jao, ao,jao) c c----------------------------------------------------------------------- c c on entry : c--------- c a, ja, ia = matrix in csr format. note that the c algorithm is in place: ao, jao can be the same c as a, ja, in which case it will be overwritten on it c upon return. c c on return : c----------- c c ao, jao = sparse matrix in modified sparse row storage format: c + ao(1:n) contains the diagonal of the matrix. c + ao(n+2:nnz) contains the nondiagonal elements of the c matrix, stored rowwise. c + jao(n+2:nnz) : their column indices c + jao(1:n+1) contains the pointer array for the nondiagonal c elements in ao(n+1:nnz) and jao(n+2:nnz). c i.e., for i .le. n+1 jao(i) points to beginning of row i c in arrays ao, jao. c here nnz = number of nonzero elements+1 c work arrays: c------------ c wk = real work array of length n c iwk = integer work array of length n+1 c c notes: c------- c Algorithm is in place. i.e. both: c c call csrmsr (n, a, ja, ia, ao, jao, ao,jao) c (in which ao, jao, are different from a, ja) c and c call csrmsr (n, a, ja, ia, a, ja, wk,iwk) c (in which wk, jwk, are different from a, ja) c are OK. c-------- c coded by Y. Saad Sep. 1989. Rechecked Feb 27, 1990. c----------------------------------------------------------------------- icount = 0 c c store away diagonal elements and count nonzero diagonal elements. c do 1 i=1,n wk(i) = 0.0d0 iwk(i+1) = ia(i+1)-ia(i) do 2 k=ia(i),ia(i+1)-1 if (ja(k) .eq. i) then wk(i) = a(k) icount = icount + 1 iwk(i+1) = iwk(i+1)-1 endif 2 continue 1 continue c c compute total length c iptr = n + ia(n+1) - icount c c copy backwards (to avoid collisions) c do 500 ii=n,1,-1 do 100 k=ia(ii+1)-1,ia(ii),-1 j = ja(k) if (j .ne. ii) then ao(iptr) = a(k) jao(iptr) = j iptr = iptr-1 endif 100 continue 500 continue c c compute pointer values and copy wk(*) c jao(1) = n+2 do 600 i=1,n ao(i) = wk(i) jao(i+1) = jao(i)+iwk(i+1) 600 continue return c------------ end of subroutine csrmsr --------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine msrcsr (n,a,ja,ao,jao,iao,wk,iwk) real*8 a(*),ao(*),wk(n) integer ja(*),jao(*),iao(n+1),iwk(n+1) c----------------------------------------------------------------------- c Modified - Sparse Row to Compressed Sparse Row c c----------------------------------------------------------------------- c converts a compressed matrix using a separated diagonal c (modified sparse row format) in the Compressed Sparse Row c format. c does not check for zero elements in the diagonal. c c c on entry : c--------- c n = row dimension of matrix c a, ja = sparse matrix in msr sparse storage format c see routine csrmsr for details on data structure c c on return : c----------- c c ao,jao,iao = output matrix in csr format. c c work arrays: c------------ c wk = real work array of length n c iwk = integer work array of length n+1 c c notes: c The original version of this was NOT in place, but has c been modified by adding the vector iwk to be in place. c The original version had ja instead of iwk everywhere in c loop 500. Modified Sun 29 May 1994 by R. Bramley (Indiana). c c----------------------------------------------------------------------- logical added do 1 i=1,n wk(i) = a(i) iwk(i) = ja(i) 1 continue iwk(n+1) = ja(n+1) iao(1) = 1 iptr = 1 c--------- do 500 ii=1,n added = .false. idiag = iptr + (iwk(ii+1)-iwk(ii)) do 100 k=iwk(ii),iwk(ii+1)-1 j = ja(k) if (j .lt. ii) then ao(iptr) = a(k) jao(iptr) = j iptr = iptr+1 elseif (added) then ao(iptr) = a(k) jao(iptr) = j iptr = iptr+1 else c add diag element - only reserve a position for it. idiag = iptr iptr = iptr+1 added = .true. c then other element ao(iptr) = a(k) jao(iptr) = j iptr = iptr+1 endif 100 continue ao(idiag) = wk(ii) jao(idiag) = ii if (.not. added) iptr = iptr+1 iao(ii+1) = iptr 500 continue return c------------ end of subroutine msrcsr --------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrcsc (n,job,ipos,a,ja,ia,ao,jao,iao) integer ia(n+1),iao(n+1),ja(*),jao(*) real*8 a(*),ao(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Compressed Sparse Column c c (transposition operation) Not in place. c----------------------------------------------------------------------- c -- not in place -- c this subroutine transposes a matrix stored in a, ja, ia format. c --------------- c on entry: c---------- c n = dimension of A. c job = integer to indicate whether to fill the values (job.eq.1) of the c matrix ao or only the pattern., i.e.,ia, and ja (job .ne.1) c c ipos = starting position in ao, jao of the transposed matrix. c the iao array takes this into account (thus iao(1) is set to ipos.) c Note: this may be useful if one needs to append the data structure c of the transpose to that of A. In this case use for example c call csrcsc (n,1,ia(n+1),a,ja,ia,a,ja,ia(n+2)) c for any other normal usage, enter ipos=1. c a = real array of length nnz (nnz=number of nonzero elements in input c matrix) containing the nonzero elements. c ja = integer array of length nnz containing the column positions c of the corresponding elements in a. c ia = integer of size n+1. ia(k) contains the position in a, ja of c the beginning of the k-th row. c c on return: c ---------- c output arguments: c ao = real array of size nzz containing the "a" part of the transpose c jao = integer array of size nnz containing the column indices. c iao = integer array of size n+1 containing the "ia" index array of c the transpose. c c----------------------------------------------------------------------- call csrcsc2 (n,n,job,ipos,a,ja,ia,ao,jao,iao) end c----------------------------------------------------------------------- subroutine csrcsc2 (n,n2,job,ipos,a,ja,ia,ao,jao,iao) integer ia(n+1),iao(n2+1),ja(*),jao(*) real*8 a(*),ao(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Compressed Sparse Column c c (transposition operation) Not in place. c----------------------------------------------------------------------- c Rectangular version. n is number of rows of CSR matrix, c n2 (input) is number of columns of CSC matrix. c----------------------------------------------------------------------- c -- not in place -- c this subroutine transposes a matrix stored in a, ja, ia format. c --------------- c on entry: c---------- c n = number of rows of CSR matrix. c n2 = number of columns of CSC matrix. c job = integer to indicate whether to fill the values (job.eq.1) of the c matrix ao or only the pattern., i.e.,ia, and ja (job .ne.1) c c ipos = starting position in ao, jao of the transposed matrix. c the iao array takes this into account (thus iao(1) is set to ipos.) c Note: this may be useful if one needs to append the data structure c of the transpose to that of A. In this case use for example c call csrcsc2 (n,n,1,ia(n+1),a,ja,ia,a,ja,ia(n+2)) c for any other normal usage, enter ipos=1. c a = real array of length nnz (nnz=number of nonzero elements in input c matrix) containing the nonzero elements. c ja = integer array of length nnz containing the column positions c of the corresponding elements in a. c ia = integer of size n+1. ia(k) contains the position in a, ja of c the beginning of the k-th row. c c on return: c ---------- c output arguments: c ao = real array of size nzz containing the "a" part of the transpose c jao = integer array of size nnz containing the column indices. c iao = integer array of size n+1 containing the "ia" index array of c the transpose. c c----------------------------------------------------------------------- c----------------- compute lengths of rows of transp(A) ---------------- do 1 i=1,n2+1 iao(i) = 0 1 continue do 3 i=1, n do 2 k=ia(i), ia(i+1)-1 j = ja(k)+1 iao(j) = iao(j)+1 2 continue 3 continue c---------- compute pointers from lengths ------------------------------ iao(1) = ipos do 4 i=1,n2 iao(i+1) = iao(i) + iao(i+1) 4 continue c--------------- now do the actual copying ----------------------------- do 6 i=1,n do 62 k=ia(i),ia(i+1)-1 j = ja(k) next = iao(j) if (job .eq. 1) ao(next) = a(k) jao(next) = i iao(j) = next+1 62 continue 6 continue c-------------------------- reshift iao and leave ---------------------- do 7 i=n2,1,-1 iao(i+1) = iao(i) 7 continue iao(1) = ipos c--------------- end of csrcsc2 ---------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrlnk (n,a,ja,ia,link) real*8 a(*) integer n, ja(*), ia(n+1), link(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Linked storage format. c----------------------------------------------------------------------- c this subroutine translates a matrix stored in compressed sparse c row into one with a linked list storage format. Only the link c array needs to be obtained since the arrays a, ja, and ia may c be unchanged and carry the same meaning for the output matrix. c in other words a, ja, ia, link is the output linked list data c structure with a, ja, unchanged from input, and ia possibly c altered (in case therea re null rows in matrix). Details on c the output array link are given below. c----------------------------------------------------------------------- c Coded by Y. Saad, Feb 21, 1991. c----------------------------------------------------------------------- c c on entry: c---------- c n = integer equal to the dimension of A. c c a = real array of size nna containing the nonzero elements c ja = integer array of size nnz containing the column positions c of the corresponding elements in a. c ia = integer of size n+1 containing the pointers to the beginning c of each row. ia(k) contains the position in a, ja of the c beginning of the k-th row. c c on return: c---------- c a, ja, are not changed. c ia may be changed if there are null rows. c c a = nonzero elements. c ja = column positions. c ia = ia(i) points to the first element of row i in linked structure. c link = integer array of size containing the linked list information. c link(k) points to the next element of the row after element c a(k), ja(k). if link(k) = 0, then there is no next element, c i.e., a(k), jcol(k) is the last element of the current row. c c Thus row number i can be accessed as follows: c next = ia(i) c while(next .ne. 0) do c value = a(next) ! value a(i,j) c jcol = ja(next) ! column index j c next = link(next) ! address of next element in row c endwhile c notes: c ------ ia may be altered on return. c----------------------------------------------------------------------- c local variables integer i, k c c loop through all rows c do 100 i =1, n istart = ia(i) iend = ia(i+1)-1 if (iend .gt. istart) then do 99 k=istart, iend-1 link(k) = k+1 99 continue link(iend) = 0 else ia(i) = 0 endif 100 continue c return c-------------end-of-csrlnk -------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine lnkcsr (n, a, jcol, istart, link, ao, jao, iao) real*8 a(*), ao(*) integer n, jcol(*), istart(n), link(*), jao(*), iao(*) c----------------------------------------------------------------------- c Linked list storage format to Compressed Sparse Row format c----------------------------------------------------------------------- c this subroutine translates a matrix stored in linked list storage c format into the compressed sparse row format. c----------------------------------------------------------------------- c Coded by Y. Saad, Feb 21, 1991. c----------------------------------------------------------------------- c c on entry: c---------- c n = integer equal to the dimension of A. c c a = real array of size nna containing the nonzero elements c jcol = integer array of size nnz containing the column positions c of the corresponding elements in a. c istart= integer array of size n poiting to the beginning of the rows. c istart(i) contains the position of the first element of c row i in data structure. (a, jcol, link). c if a row is empty istart(i) must be zero. c link = integer array of size nnz containing the links in the linked c list data structure. link(k) points to the next element c of the row after element ao(k), jcol(k). if link(k) = 0, c then there is no next element, i.e., ao(k), jcol(k) is c the last element of the current row. c c on return: c----------- c ao, jao, iao = matrix stored in csr format: c c ao = real array containing the values of the nonzero elements of c the matrix stored row-wise. c jao = integer array of size nnz containing the column indices. c iao = integer array of size n+1 containing the pointers array to the c beginning of each row. iao(i) is the address in ao,jao of c first element of row i. c c----------------------------------------------------------------------- c first determine individial bandwidths and pointers. c----------------------------------------------------------------------- c local variables integer irow, ipos, next c----------------------------------------------------------------------- ipos = 1 iao(1) = ipos c c loop through all rows c do 100 irow =1, n c c unroll i-th row. c next = istart(irow) 10 if (next .eq. 0) goto 99 jao(ipos) = jcol(next) ao(ipos) = a(next) ipos = ipos+1 next = link(next) goto 10 99 iao(irow+1) = ipos 100 continue c return c-------------end-of-lnkcsr ------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrdia (n,idiag,job,a,ja,ia,ndiag, * diag,ioff,ao,jao,iao,ind) real*8 diag(ndiag,idiag), a(*), ao(*) integer ia(*), ind(*), ja(*), jao(*), iao(*), ioff(*) c----------------------------------------------------------------------- c Compressed sparse row to diagonal format c----------------------------------------------------------------------- c this subroutine extracts idiag diagonals from the input matrix a, c a, ia, and puts the rest of the matrix in the output matrix ao, c jao, iao. The diagonals to be extracted depend on the value of job c (see below for details.) In the first case, the diagonals to be c extracted are simply identified by their offsets provided in ioff c by the caller. In the second case, the code internally determines c the idiag most significant diagonals, i.e., those diagonals of the c matrix which have the largest number of nonzero elements, and c extracts them. c----------------------------------------------------------------------- c on entry: c---------- c n = dimension of the matrix a. c idiag = integer equal to the number of diagonals to be extracted. c Note: on return idiag may be modified. c a, ja, c ia = matrix stored in a, ja, ia, format c job = integer. serves as a job indicator. Job is better thought c of as a two-digit number job=xy. If the first (x) digit c is one on entry then the diagonals to be extracted are c internally determined. In this case csrdia exctracts the c idiag most important diagonals, i.e. those having the largest c number on nonzero elements. If the first digit is zero c then csrdia assumes that ioff(*) contains the offsets c of the diagonals to be extracted. there is no verification c that ioff(*) contains valid entries. c The second (y) digit of job determines whether or not c the remainder of the matrix is to be written on ao,jao,iao. c If it is zero then ao, jao, iao is not filled, i.e., c the diagonals are found and put in array diag and the rest is c is discarded. if it is one, ao, jao, iao contains matrix c of the remaining elements. c Thus: c job= 0 means do not select diagonals internally (pick those c defined by ioff) and do not fill ao,jao,iao c job= 1 means do not select diagonals internally c and fill ao,jao,iao c job=10 means select diagonals internally c and do not fill ao,jao,iao c job=11 means select diagonals internally c and fill ao,jao,iao c c ndiag = integer equal to the first dimension of array diag. c c on return: c----------- c c idiag = number of diagonals found. This may be smaller than its value c on entry. c diag = real array of size (ndiag x idiag) containing the diagonals c of A on return c c ioff = integer array of length idiag, containing the offsets of the c diagonals to be extracted. c ao, jao c iao = remainder of the matrix in a, ja, ia format. c work arrays: c------------ c ind = integer array of length 2*n-1 used as integer work space. c needed only when job.ge.10 i.e., in case the diagonals are to c be selected internally. c c Notes: c------- c 1) The algorithm is in place: ao, jao, iao can be overwritten on c a, ja, ia if desired c 2) When the code is required to select the diagonals (job .ge. 10) c the selection of the diagonals is done from left to right c as a result if several diagonals have the same weight (number c of nonzero elemnts) the leftmost one is selected first. c----------------------------------------------------------------------- job1 = job/10 job2 = job-job1*10 if (job1 .eq. 0) goto 50 n2 = n+n-1 call infdia(n,ja,ia,ind,idum) c----------- determine diagonals to accept.---------------------------- c----------------------------------------------------------------------- ii = 0 4 ii=ii+1 jmax = 0 do 41 k=1, n2 j = ind(k) if (j .le. jmax) goto 41 i = k jmax = j 41 continue if (jmax .le. 0) then ii = ii-1 goto 42 endif ioff(ii) = i-n ind(i) = - jmax if (ii .lt. idiag) goto 4 42 idiag = ii c---------------- initialize diago to zero ----------------------------- 50 continue do 55 j=1,idiag do 54 i=1,n diag(i,j) = 0.0d0 54 continue 55 continue c----------------------------------------------------------------------- ko = 1 c----------------------------------------------------------------------- c extract diagonals and accumulate remaining matrix. c----------------------------------------------------------------------- do 6 i=1, n do 51 k=ia(i),ia(i+1)-1 j = ja(k) do 52 l=1,idiag if (j-i .ne. ioff(l)) goto 52 diag(i,l) = a(k) goto 51 52 continue c--------------- append element not in any diagonal to ao,jao,iao ----- if (job2 .eq. 0) goto 51 ao(ko) = a(k) jao(ko) = j ko = ko+1 51 continue if (job2 .ne. 0 ) ind(i+1) = ko 6 continue if (job2 .eq. 0) return c finish with iao iao(1) = 1 do 7 i=2,n+1 iao(i) = ind(i) 7 continue return c----------- end of csrdia --------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine diacsr (n,job,idiag,diag,ndiag,ioff,a,ja,ia) real*8 diag(ndiag,idiag), a(*), t integer ia(*), ja(*), ioff(*) c----------------------------------------------------------------------- c diagonal format to compressed sparse row c----------------------------------------------------------------------- c this subroutine extract the idiag most important diagonals from the c input matrix a, ja, ia, i.e, those diagonals of the matrix which have c the largest number of nonzero elements. If requested (see job), c the rest of the matrix is put in a the output matrix ao, jao, iao c----------------------------------------------------------------------- c on entry: c---------- c n = integer. dimension of the matrix a. c job = integer. job indicator with the following meaning. c if (job .eq. 0) then check for each entry in diag c whether this entry is zero. If it is then do not include c in the output matrix. Note that the test is a test for c an exact arithmetic zero. Be sure that the zeros are c actual zeros in double precision otherwise this would not c work. c c idiag = integer equal to the number of diagonals to be extracted. c Note: on return idiag may be modified. c c diag = real array of size (ndiag x idiag) containing the diagonals c of A on return. c c ndiag = integer equal to the first dimension of array diag. c c ioff = integer array of length idiag, containing the offsets of the c diagonals to be extracted. c c on return: c----------- c a, c ja, c ia = matrix stored in a, ja, ia, format c c Note: c ----- the arrays a and ja should be of length n*idiag. c c----------------------------------------------------------------------- ia(1) = 1 ko = 1 do 80 i=1, n do 70 jj = 1, idiag j = i+ioff(jj) if (j .lt. 1 .or. j .gt. n) goto 70 t = diag(i,jj) if (job .eq. 0 .and. t .eq. 0.0d0) goto 70 a(ko) = t ja(ko) = j ko = ko+1 70 continue ia(i+1) = ko 80 continue return c----------- end of diacsr --------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine bsrcsr (job, n, m, na, a, ja, ia, ao, jao, iao) implicit none integer job, n, m, na, ia(*), ja(*), jao(*), iao(n+1) real*8 a(na,*), ao(*) c----------------------------------------------------------------------- c Block Sparse Row to Compressed Sparse Row. c----------------------------------------------------------------------- c NOTE: ** meanings of parameters may have changed wrt earlier versions c FORMAT DEFINITION HAS CHANGED WRT TO EARLIER VERSIONS... c----------------------------------------------------------------------- c c converts a matrix stored in block-reduced a, ja, ia format to the c general sparse row a, ja, ia format. A matrix that has a block c structure is a matrix whose entries are blocks of the same size m c (e.g. 3 x 3). Then it is often preferred to work with the reduced c graph of the matrix. Instead of storing one element at a time one can c store a whole block at a time. In this storage scheme an entry is a c square array holding the m**2 elements of a block. c c----------------------------------------------------------------------- c on entry: c---------- c job = if job.eq.0 on entry, values are not copied (pattern only) c c n = the block row dimension of the matrix. c c m = the dimension of each block. Thus, the actual row dimension c of A is n x m. c c na = first dimension of array a as declared in calling program. c This should be .ge. m**2. c c a = real array containing the real entries of the matrix. Recall c that each entry is in fact an m x m block. These entries c are stored column-wise in locations a(1:m*m,k) for each k-th c entry. See details below. c c ja = integer array of length n. ja(k) contains the column index c of the leading element, i.e., the element (1,1) of the block c that is held in the column a(*,k) of the value array. c c ia = integer array of length n+1. ia(i) points to the beginning c of block row number i in the arrays a and ja. c c on return: c----------- c ao, jao, c iao = matrix stored in compressed sparse row format. The number of c rows in the new matrix is n x m. c c Notes: THIS CODE IS NOT IN PLACE. c c----------------------------------------------------------------------- c BSR FORMAT. c---------- c Each row of A contains the m x m block matrix unpacked column- c wise (this allows the user to declare the array a as a(m,m,*) on entry c if desired). The block rows are stored in sequence just as for the c compressed sparse row format. c c----------------------------------------------------------------------- c example with m = 2: c 1 2 3 c +-------|--------|--------+ +-------+ c | 1 2 | 0 0 | 3 4 | Block | x 0 x | 1 c | 5 6 | 0 0 | 7 8 | Representation: | 0 x x | 2 c +-------+--------+--------+ | x 0 0 | 3 c | 0 0 | 9 10 | 11 12 | +-------+ c | 0 0 | 13 14 | 15 16 | c +-------+--------+--------+ c | 17 18 | 0 0 | 0 0 | c | 22 23 | 0 0 | 0 0 | c +-------+--------+--------+ c c For this matrix: n = 3 c m = 2 c nnz = 5 c----------------------------------------------------------------------- c Data structure in Block Sparse Row format: c------------------------------------------- c Array A: c------------------------- c 1 3 9 11 17 <<--each m x m block is stored column-wise c 5 7 13 15 22 in a column of the array A. c 2 4 10 12 18 c 6 8 14 16 23 c------------------------- c JA 1 3 2 3 1 <<-- column indices for each block. Note that c------------------------- these indices are wrt block matrix. c IA 1 3 5 6 <<-- pointers to beginning of each block row c------------------------- in arrays A and JA. c----------------------------------------------------------------------- c locals c integer i, i1, i2, ij, ii, irow, j, jstart, k, krow, no logical val c val = (job.ne.0) no = n * m irow = 1 krow = 1 iao(irow) = 1 c----------------------------------------------------------------------- do 2 ii=1, n c c recall: n is the block-row dimension c i1 = ia(ii) i2 = ia(ii+1)-1 c c create m rows for each block row -- i.e., each k. c do 23 i=1,m do 21 k=i1, i2 jstart = m*(ja(k)-1) do 22 j=1,m ij = (j-1)*m + i if (val) ao(krow) = a(ij,k) jao(krow) = jstart+j krow = krow+1 22 continue 21 continue irow = irow+1 iao(irow) = krow 23 continue 2 continue return c-------------end-of-bsrcsr -------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrbsr (job,nrow,m,na,a,ja,ia,ao,jao,iao,iw,ierr) implicit none integer job,ierr,nrow,m,na,ia(nrow+1),ja(*),jao(na),iao(*),iw(*) real*8 a(*),ao(na,*) c----------------------------------------------------------------------- c Compressed Sparse Row to Block Sparse Row c----------------------------------------------------------------------- c c This subroutine converts a matrix stored in a general compressed a, c ja, ia format into a a block sparse row format a(m,m,*),ja(*),ia(*). c See routine bsrcsr for more details on data structure for block c matrices. c c NOTES: 1) the initial matrix does not have to have a block structure. c zero padding is done for general sparse matrices. c 2) For most practical purposes, na should be the same as m*m. c c----------------------------------------------------------------------- c c In what follows nr=1+(nrow-1)/m = block-row dimension of output matrix c c on entry: c---------- c c job = job indicator. c job = 0 -> only the pattern of output matrix is generated c job > 0 -> both pattern and values are generated. c job = -1 -> iao(1) will return the number of nonzero blocks, c in the output matrix. In this case jao(1:nr) is used as c workspace, ao is untouched, iao is untouched except iao(1) c c nrow = integer, the actual row dimension of the matrix. c c m = integer equal to the dimension of each block. m should be > 0. c c na = first dimension of array ao as declared in calling program. c na should be .ge. m*m. c c a, ja, c ia = input matrix stored in compressed sparse row format. c c on return: c----------- c c ao = real array containing the values of the matrix. For details c on the format see below. Each row of a contains the m x m c block matrix unpacked column-wise (this allows the user to c declare the array a as ao(m,m,*) on entry if desired). The c block rows are stored in sequence just as for the compressed c sparse row format. The block dimension of the output matrix c is nr = 1 + (nrow-1) / m. c c jao = integer array. containing the block-column indices of the c block-matrix. Each jao(k) is an integer between 1 and nr c containing the block column index of the block ao(*,k). c c iao = integer array of length nr+1. iao(i) points to the beginning c of block row number i in the arrays ao and jao. When job=-1 c iao(1) contains the number of nonzero blocks of the output c matrix and the rest of iao is unused. This is useful for c determining the lengths of ao and jao. c c ierr = integer, error code. c 0 -- normal termination c 1 -- m is equal to zero c 2 -- NA too small to hold the blocks (should be .ge. m**2) c c Work arrays: c------------- c iw = integer work array of dimension nr = 1 + (nrow-1) / m c c NOTES: c------- c 1) this code is not in place. c 2) see routine bsrcsr for details on data sctructure for block c sparse row format. c c----------------------------------------------------------------------- c nr is the block-dimension of the output matrix. c integer nr, m2, io, ko, ii, len, k, jpos, j, i, ij, jr, irow logical vals c----- ierr = 0 if (m*m .gt. na) ierr = 2 if (m .eq. 0) ierr = 1 if (ierr .ne. 0) return c----------------------------------------------------------------------- vals = (job .gt. 0) nr = 1 + (nrow-1) / m m2 = m*m ko = 1 io = 1 iao(io) = 1 len = 0 c c iw determines structure of block-row (nonzero indicator) c do j=1, nr iw(j) = 0 enddo c c big loop -- leap by m rows each time. c do ii=1, nrow, m irow = 0 c c go through next m rows -- make sure not to go beyond nrow. c do while (ii+irow .le. nrow .and. irow .le. m-1) do k=ia(ii+irow),ia(ii+irow+1)-1 c c block column index = (scalar column index -1) / m + 1 c j = ja(k)-1 jr = j/m + 1 j = j - (jr-1)*m jpos = iw(jr) if (jpos .eq. 0) then c c create a new block c iw(jr) = ko jao(ko) = jr if (vals) then c c initialize new block to zero -- then copy nonzero element c do i=1, m2 ao(i,ko) = 0.0d0 enddo ij = j*m + irow + 1 ao(ij,ko) = a(k) endif ko = ko+1 else c c copy column index and nonzero element c jao(jpos) = jr ij = j*m + irow + 1 if (vals) ao(ij,jpos) = a(k) endif enddo irow = irow+1 enddo c c refresh iw c do j = iao(io),ko-1 iw(jao(j)) = 0 enddo if (job .eq. -1) then len = len + ko-1 ko = 1 else io = io+1 iao(io) = ko endif enddo if (job .eq. -1) iao(1) = len c return c--------------end-of-csrbsr-------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrbnd (n,a,ja,ia,job,abd,nabd,lowd,ml,mu,ierr) real*8 a(*),abd(nabd,n) integer ia(n+1),ja(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Banded (Linpack ) format. c----------------------------------------------------------------------- c this subroutine converts a general sparse matrix stored in c compressed sparse row format into the banded format. for the c banded format,the Linpack conventions are assumed (see below). c----------------------------------------------------------------------- c on entry: c---------- c n = integer,the actual row dimension of the matrix. c c a, c ja, c ia = input matrix stored in compressed sparse row format. c c job = integer. if job=1 then the values of the lower bandwith ml c and the upper bandwidth mu are determined internally. c otherwise it is assumed that the values of ml and mu c are the correct bandwidths on input. See ml and mu below. c c nabd = integer. first dimension of array abd. c c lowd = integer. this should be set to the row number in abd where c the lowest diagonal (leftmost) of A is located. c lowd should be ( 1 .le. lowd .le. nabd). c if it is not known in advance what lowd should be c enter lowd = 0 and the default value lowd = ml+mu+1 c will be chosen. Alternative: call routine getbwd from unary c first to detrermione ml and mu then define lowd accordingly. c (Note: the banded solvers in linpack use lowd=2*ml+mu+1. ) c c ml = integer. equal to the bandwidth of the strict lower part of A c mu = integer. equal to the bandwidth of the strict upper part of A c thus the total bandwidth of A is ml+mu+1. c if ml+mu+1 is found to be larger than lowd then an error c flag is raised (unless lowd = 0). see ierr. c c note: ml and mu are assumed to have the correct bandwidth values c as defined above if job is set to zero on entry. c c on return: c----------- c c abd = real array of dimension abd(nabd,n). c on return contains the values of the matrix stored in c banded form. The j-th column of abd contains the elements c of the j-th column of the original matrix comprised in the c band ( i in (j-ml,j+mu) ) with the lowest diagonal at c the bottom row (row lowd). See details below for this format. c c ml = integer. equal to the bandwidth of the strict lower part of A c mu = integer. equal to the bandwidth of the strict upper part of A c if job=1 on entry then these two values are internally computed. c c lowd = integer. row number in abd where the lowest diagonal c (leftmost) of A is located on return. In case lowd = 0 c on return, then it is defined to ml+mu+1 on return and the c lowd will contain this value on return. ` c c ierr = integer. used for error messages. On return: c ierr .eq. 0 :means normal return c ierr .eq. -1 : means invalid value for lowd. (either .lt. 0 c or larger than nabd). c ierr .eq. -2 : means that lowd is not large enough and as c result the matrix cannot be stored in array abd. c lowd should be at least ml+mu+1, where ml and mu are as c provided on output. c c----------------------------------------------------------------------* c Additional details on banded format. (this closely follows the * c format used in linpack. may be useful for converting a matrix into * c this storage format in order to use the linpack banded solvers). * c----------------------------------------------------------------------* c --- band storage format for matrix abd --- * c uses ml+mu+1 rows of abd(nabd,*) to store the diagonals of * c a in rows of abd starting from the lowest (sub)-diagonal which is * c stored in row number lowd of abd. the minimum number of rows needed * c in abd is ml+mu+1, i.e., the minimum value for lowd is ml+mu+1. the * c j-th column of abd contains the elements of the j-th column of a, * c from bottom to top: the element a(j+ml,j) is stored in position * c abd(lowd,j), then a(j+ml-1,j) in position abd(lowd-1,j) and so on. * c Generally, the element a(j+k,j) of original matrix a is stored in * c position abd(lowd+k-ml,j), for k=ml,ml-1,..,0,-1, -mu. * c The first dimension nabd of abd must be .ge. lowd * c * c example [from linpack ]: if the original matrix is * c * c 11 12 13 0 0 0 * c 21 22 23 24 0 0 * c 0 32 33 34 35 0 original banded matrix * c 0 0 43 44 45 46 * c 0 0 0 54 55 56 * c 0 0 0 0 65 66 * c * c then n = 6, ml = 1, mu = 2. lowd should be .ge. 4 (=ml+mu+1) and * c if lowd = 5 for example, abd should be: * c * c untouched --> x x x x x x * c * * 13 24 35 46 * c * 12 23 34 45 56 resulting abd matrix in banded * c 11 22 33 44 55 66 format * c row lowd--> 21 32 43 54 65 * * c * c * = not used * c * c----------------------------------------------------------------------* c first determine ml and mu. c----------------------------------------------------------------------- ierr = 0 c----------- if (job .eq. 1) call getbwd(n,a,ja,ia,ml,mu) m = ml+mu+1 if (lowd .eq. 0) lowd = m if (m .gt. lowd) ierr = -2 if (lowd .gt. nabd .or. lowd .lt. 0) ierr = -1 if (ierr .lt. 0) return c------------ do 15 i=1,m ii = lowd -i+1 do 10 j=1,n abd(ii,j) = 0.0d0 10 continue 15 continue c--------------------------------------------------------------------- mdiag = lowd-ml do 30 i=1,n do 20 k=ia(i),ia(i+1)-1 j = ja(k) abd(i-j+mdiag,j) = a(k) 20 continue 30 continue return c------------- end of csrbnd ------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine bndcsr (n,abd,nabd,lowd,ml,mu,a,ja,ia,len,ierr) real*8 a(*),abd(nabd,*), t integer ia(n+1),ja(*) c----------------------------------------------------------------------- c Banded (Linpack ) format to Compressed Sparse Row format. c----------------------------------------------------------------------- c on entry: c---------- c n = integer,the actual row dimension of the matrix. c c nabd = first dimension of array abd. c c abd = real array containing the values of the matrix stored in c banded form. The j-th column of abd contains the elements c of the j-th column of the original matrix,comprised in the c band ( i in (j-ml,j+mu) ) with the lowest diagonal located c in row lowd (see below). c c lowd = integer. this should be set to the row number in abd where c the lowest diagonal (leftmost) of A is located. c lowd should be s.t. ( 1 .le. lowd .le. nabd). c The subroutines dgbco, ... of linpack use lowd=2*ml+mu+1. c c ml = integer. equal to the bandwidth of the strict lower part of A c mu = integer. equal to the bandwidth of the strict upper part of A c thus the total bandwidth of A is ml+mu+1. c if ml+mu+1 is found to be larger than nabd then an error c message is set. see ierr. c c len = integer. length of arrays a and ja. bndcsr will stop if the c length of the arrays a and ja is insufficient to store the c matrix. see ierr. c c on return: c----------- c a, c ja, c ia = input matrix stored in compressed sparse row format. c c lowd = if on entry lowd was zero then lowd is reset to the default c value ml+mu+l. c c ierr = integer. used for error message output. c ierr .eq. 0 :means normal return c ierr .eq. -1 : means invalid value for lowd. c ierr .gt. 0 : means that there was not enough storage in a and ja c for storing the ourput matrix. The process ran out of space c (as indicated by len) while trying to fill row number ierr. c This should give an idea of much more storage might be required. c Moreover, the first irow-1 rows are correctly filled. c c notes: the values in abd found to be equal to zero c ----- (actual test: if (abd(...) .eq. 0.0d0) are removed. c The resulting may not be identical to a csr matrix c originally transformed to a bnd format. c c----------------------------------------------------------------------- ierr = 0 c----------- if (lowd .gt. nabd .or. lowd .le. 0) then ierr = -1 return endif c----------- ko = 1 ia(1) = 1 do 30 irow=1,n c----------------------------------------------------------------------- i = lowd do 20 j=irow-ml,irow+mu if (j .le. 0 ) goto 19 if (j .gt. n) goto 21 t = abd(i,j) if (t .eq. 0.0d0) goto 19 if (ko .gt. len) then ierr = irow return endif a(ko) = t ja(ko) = j ko = ko+1 19 i = i-1 20 continue c end for row irow 21 ia(irow+1) = ko 30 continue return c------------- end of bndcsr ------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrssk (n,imod,a,ja,ia,asky,isky,nzmax,ierr) real*8 a(*),asky(nzmax) integer n, imod, nzmax, ierr, ia(n+1), isky(n+1), ja(*) c----------------------------------------------------------------------- c Compressed Sparse Row to Symmetric Skyline Format c or Symmetric Sparse Row c----------------------------------------------------------------------- c this subroutine translates a compressed sparse row or a symmetric c sparse row format into a symmetric skyline format. c the input matrix can be in either compressed sparse row or the c symmetric sparse row format. The output matrix is in a symmetric c skyline format: a real array containing the (active portions) of the c rows in sequence and a pointer to the beginning of each row. c c This module is NOT in place. c----------------------------------------------------------------------- c Coded by Y. Saad, Oct 5, 1989. Revised Feb. 18, 1991. c----------------------------------------------------------------------- c c on entry: c---------- c n = integer equal to the dimension of A. c imod = integer indicating the variant of skyline format wanted: c imod = 0 means the pointer isky points to the `zeroth' c element of the row, i.e., to the position of the diagonal c element of previous row (for i=1, isky(1)= 0) c imod = 1 means that itpr points to the beginning of the row. c imod = 2 means that isky points to the end of the row (diagonal c element) c c a = real array of size nna containing the nonzero elements c ja = integer array of size nnz containing the column positions c of the corresponding elements in a. c ia = integer of size n+1. ia(k) contains the position in a, ja of c the beginning of the k-th row. c nzmax = integer. must be set to the number of available locations c in the output array asky. c c on return: c---------- c c asky = real array containing the values of the matrix stored in skyline c format. asky contains the sequence of active rows from c i=1, to n, an active row being the row of elemnts of c the matrix contained between the leftmost nonzero element c and the diagonal element. c isky = integer array of size n+1 containing the pointer array to c each row. The meaning of isky depends on the input value of c imod (see above). c ierr = integer. Error message. If the length of the c output array asky exceeds nzmax. ierr returns the minimum value c needed for nzmax. otherwise ierr=0 (normal return). c c Notes: c 1) This module is NOT in place. c 2) even when imod = 2, length of isky is n+1, not n. c c----------------------------------------------------------------------- c first determine individial bandwidths and pointers. c----------------------------------------------------------------------- ierr = 0 isky(1) = 0 do 3 i=1,n ml = 0 do 31 k=ia(i),ia(i+1)-1 ml = max(ml,i-ja(k)+1) 31 continue isky(i+1) = isky(i)+ml 3 continue c c test if there is enough space asky to do the copying. c nnz = isky(n+1) if (nnz .gt. nzmax) then ierr = nnz return endif c c fill asky with zeros. c do 1 k=1, nnz asky(k) = 0.0d0 1 continue c c copy nonzero elements. c do 4 i=1,n kend = isky(i+1) do 41 k=ia(i),ia(i+1)-1 j = ja(k) if (j .le. i) asky(kend+j-i) = a(k) 41 continue 4 continue c c modify pointer according to imod if necessary. c if (imod .eq. 0) return if (imod .eq. 1) then do 50 k=1, n+1 isky(k) = isky(k)+1 50 continue endif if (imod .eq. 2) then do 60 k=1, n isky(k) = isky(k+1) 60 continue endif c return c------------- end of csrssk ------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine sskssr (n,imod,asky,isky,ao,jao,iao,nzmax,ierr) real*8 asky(*),ao(nzmax) integer n, imod,nzmax,ierr, isky(n+1),iao(n+1),jao(nzmax) c----------------------------------------------------------------------- c Symmetric Skyline Format to Symmetric Sparse Row format. c----------------------------------------------------------------------- c tests for exact zeros in skyline matrix (and ignores them in c output matrix). In place routine (a, isky :: ao, iao) c----------------------------------------------------------------------- c this subroutine translates a symmetric skyline format into a c symmetric sparse row format. Each element is tested to see if it is c a zero element. Only the actual nonzero elements are retained. Note c that the test used is simple and does take into account the smallness c of a value. the subroutine filter (see unary module) can be used c for this purpose. c----------------------------------------------------------------------- c Coded by Y. Saad, Oct 5, 1989. Revised Feb 18, 1991./ c----------------------------------------------------------------------- c c on entry: c---------- c n = integer equal to the dimension of A. c imod = integer indicating the variant of skyline format used: c imod = 0 means the pointer iao points to the `zeroth' c element of the row, i.e., to the position of the diagonal c element of previous row (for i=1, iao(1)= 0) c imod = 1 means that itpr points to the beginning of the row. c imod = 2 means that iao points to the end of the row c (diagonal element) c asky = real array containing the values of the matrix. asky contains c the sequence of active rows from i=1, to n, an active row c being the row of elemnts of the matrix contained between the c leftmost nonzero element and the diagonal element. c isky = integer array of size n+1 containing the pointer array to c each row. isky (k) contains the address of the beginning of the c k-th active row in the array asky. c nzmax = integer. equal to the number of available locations in the c output array ao. c c on return: c ---------- c ao = real array of size nna containing the nonzero elements c jao = integer array of size nnz containing the column positions c of the corresponding elements in a. c iao = integer of size n+1. iao(k) contains the position in a, ja of c the beginning of the k-th row. c ierr = integer. Serving as error message. If the length of the c output arrays ao, jao exceeds nzmax then ierr returns c the row number where the algorithm stopped: rows c i, to ierr-1 have been processed succesfully. c ierr = 0 means normal return. c ierr = -1 : illegal value for imod c Notes: c------- c This module is in place: ao and iao can be the same as asky, and isky. c----------------------------------------------------------------------- c local variables integer next, kend, kstart, i, j ierr = 0 c c check for validity of imod c if (imod.ne.0 .and. imod.ne.1 .and. imod .ne. 2) then ierr =-1 return endif c c next = pointer to next available position in output matrix c kend = pointer to end of current row in skyline matrix. c next = 1 c c set kend = start position -1 in skyline matrix. c kend = 0 if (imod .eq. 1) kend = isky(1)-1 if (imod .eq. 0) kend = isky(1) c c loop through all rows c do 50 i=1,n c c save value of pointer to ith row in output matrix c iao(i) = next c c get beginnning and end of skyline row c kstart = kend+1 if (imod .eq. 0) kend = isky(i+1) if (imod .eq. 1) kend = isky(i+1)-1 if (imod .eq. 2) kend = isky(i) c c copy element into output matrix unless it is a zero element. c do 40 k=kstart,kend if (asky(k) .eq. 0.0d0) goto 40 j = i-(kend-k) jao(next) = j ao(next) = asky(k) next=next+1 if (next .gt. nzmax+1) then ierr = i return endif 40 continue 50 continue iao(n+1) = next return c-------------end-of-sskssr -------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine csrjad (nrow, a, ja, ia, idiag, iperm, ao, jao, iao) integer ja(*), jao(*), ia(nrow+1), iperm(nrow), iao(nrow) real*8 a(*), ao(*) c----------------------------------------------------------------------- c Compressed Sparse Row to JAgged Diagonal storage. c----------------------------------------------------------------------- c this subroutine converts matrix stored in the compressed sparse c row format to the jagged diagonal format. The data structure c for the JAD (Jagged Diagonal storage) is as follows. The rows of c the matrix are (implicitly) permuted so that their lengths are in c decreasing order. The real entries ao(*) and their column indices c jao(*) are stored in succession. The number of such diagonals is idiag. c the lengths of each of these diagonals is stored in iao(*). c For more details see [E. Anderson and Y. Saad, c ``Solving sparse triangular systems on parallel computers'' in c Inter. J. of High Speed Computing, Vol 1, pp. 73-96 (1989).] c or [Y. Saad, ``Krylov Subspace Methods on Supercomputers'' c SIAM J. on Stat. Scient. Comput., volume 10, pp. 1200-1232 (1989).] c----------------------------------------------------------------------- c on entry: c---------- c nrow = row dimension of the matrix A. c c a, c ia, c ja = input matrix in compressed sparse row format. c c on return: c---------- c c idiag = integer. The number of jagged diagonals in the matrix. c c iperm = integer array of length nrow containing the permutation c of the rows that leads to a decreasing order of the c number of nonzero elements. c c ao = real array containing the values of the matrix A in c jagged diagonal storage. The j-diagonals are stored c in ao in sequence. c c jao = integer array containing the column indices of the c entries in ao. c c iao = integer array containing pointers to the beginning c of each j-diagonal in ao, jao. iao is also used as c a work array and it should be of length n at least. c c----------------------------------------------------------------------- c ---- define initial iperm and get lengths of each row c ---- jao is used a work vector to store tehse lengths c idiag = 0 ilo = nrow do 10 j=1, nrow iperm(j) = j len = ia(j+1) - ia(j) ilo = min(ilo,len) idiag = max(idiag,len) jao(j) = len 10 continue c c call sorter to get permutation. use iao as work array. c call dcsort (jao, nrow, iao, iperm, ilo, idiag) c c define output data structure. first lengths of j-diagonals c do 20 j=1, nrow iao(j) = 0 20 continue do 40 k=1, nrow len = jao(iperm(k)) do 30 i=1,len iao(i) = iao(i)+1 30 continue 40 continue c c get the output matrix itself c k1 = 1 k0 = k1 do 60 jj=1, idiag len = iao(jj) do 50 k=1,len i = ia(iperm(k))+jj-1 ao(k1) = a(i) jao(k1) = ja(i) k1 = k1+1 50 continue iao(jj) = k0 k0 = k1 60 continue iao(idiag+1) = k1 return c----------end-of-csrjad------------------------------------------------ c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine jadcsr (nrow, idiag, a, ja, ia, iperm, ao, jao, iao) integer ja(*), jao(*), ia(idiag+1), iperm(nrow), iao(nrow+1) real*8 a(*), ao(*) c----------------------------------------------------------------------- c Jagged Diagonal Storage to Compressed Sparse Row c----------------------------------------------------------------------- c this subroutine converts a matrix stored in the jagged diagonal format c to the compressed sparse row format. c----------------------------------------------------------------------- c on entry: c---------- c nrow = integer. the row dimension of the matrix A. c c idiag = integer. The number of jagged diagonals in the data c structure a, ja, ia. c c a, c ja, c ia = input matrix in jagged diagonal format. c c iperm = permutation of the rows used to obtain the JAD ordering. c c on return: c---------- c c ao, jao, c iao = matrix in CSR format. c----------------------------------------------------------------------- c determine first the pointers for output matrix. Go through the c structure once: c do 137 j=1,nrow jao(j) = 0 137 continue c c compute the lengths of each row of output matrix - c do 140 i=1, idiag len = ia(i+1)-ia(i) do 138 k=1,len jao(iperm(k)) = jao(iperm(k))+1 138 continue 140 continue c c remember to permute c kpos = 1 iao(1) = 1 do 141 i=1, nrow kpos = kpos+jao(i) iao(i+1) = kpos 141 continue c c copy elemnts one at a time. c do 200 jj = 1, idiag k1 = ia(jj)-1 len = ia(jj+1)-k1-1 do 160 k=1,len kpos = iao(iperm(k)) ao(kpos) = a(k1+k) jao(kpos) = ja(k1+k) iao(iperm(k)) = kpos+1 160 continue 200 continue c c rewind pointers c do 5 j=nrow,1,-1 iao(j+1) = iao(j) 5 continue iao(1) = 1 return c----------end-of-jadcsr------------------------------------------------ c----------------------------------------------------------------------- end subroutine dcsort(ival, n, icnt, index, ilo, ihi) c----------------------------------------------------------------------- c Specifications for arguments: c ---------------------------- integer n, ilo, ihi, ival(n), icnt(ilo:ihi), index(n) c----------------------------------------------------------------------- c This routine computes a permutation which, when applied to the c input vector ival, sorts the integers in ival in descending c order. The permutation is represented by the vector index. The c permuted ival can be interpreted as follows: c ival(index(i-1)) .ge. ival(index(i)) .ge. ival(index(i+1)) c c A specialized sort, the distribution counting sort, is used c which takes advantage of the knowledge that c 1) The values are in the (small) range [ ilo, ihi ] c 2) Values are likely to be repeated often c c contributed to SPARSKIT by Mike Heroux. (Cray Research) c --------------------------------------- c----------------------------------------------------------------------- c Usage: c------ c call dcsort( ival, n, icnt, index, ilo, ihi ) c c Arguments: c----------- c ival integer array (input) c On entry, ia is an n dimensional array that contains c the values to be sorted. ival is unchanged on exit. c c n integer (input) c On entry, n is the number of elements in ival and index. c c icnt integer (work) c On entry, is an integer work vector of length c (ihi - ilo + 1). c c index integer array (output) c On exit, index is an n-length integer vector containing c the permutation which sorts the vector ival. c c ilo integer (input) c On entry, ilo is .le. to the minimum value in ival. c c ihi integer (input) c On entry, ihi is .ge. to the maximum value in ival. c c Remarks: c--------- c The permutation is NOT applied to the vector ival. c c---------------------------------------------------------------- c c Local variables: c Other integer values are temporary indices. c c Author: c-------- c Michael Heroux c Sandra Carney c Mathematical Software Research Group c Cray Research, Inc. c c References: c Knuth, Donald E., "The Art of Computer Programming, Volume 3: c Sorting and Searching," Addison-Wesley, Reading, Massachusetts, c 1973, pp. 78-79. c c Revision history: c 05/09/90: Original implementation. A variation of the c Distribution Counting Sort recommended by c Sandra Carney. (Mike Heroux) c c----------------------------------------------------------------- c ---------------------------------- c Specifications for local variables c ---------------------------------- integer i, j, ivalj c c -------------------------- c First executable statement c -------------------------- do 10 i = ilo, ihi icnt(i) = 0 10 continue c do 20 i = 1, n icnt(ival(i)) = icnt(ival(i)) + 1 20 continue c do 30 i = ihi-1,ilo,-1 icnt(i) = icnt(i) + icnt(i+1) 30 continue c do 40 j = n, 1, -1 ivalj = ival(j) index(icnt(ivalj)) = j icnt(ivalj) = icnt(ivalj) - 1 40 continue return end c-------end-of-dcsort--------------------------------------------------- c----------------------------------------------------------------------- subroutine cooell(job,n,nnz,a,ja,ia,ao,jao,lda,ncmax,nc,ierr) implicit none integer job,n,nnz,lda,ncmax,nc,ierr integer ja(nnz),ia(nnz),jao(lda,ncmax) real*8 a(nnz),ao(lda,ncmax) c----------------------------------------------------------------------- c COOrdinate format to ELLpack format c----------------------------------------------------------------------- c On entry: c job -- 0 if only pattern is to be processed(AO is not touched) c n -- number of rows in the matrix c a,ja,ia -- input matix in COO format c lda -- leading dimension of array AO and JAO c ncmax -- size of the second dimension of array AO and JAO c c On exit: c ao,jao -- the matrix in ELL format c nc -- maximum number of nonzeros per row c ierr -- 0 if convertion succeeded c -1 if LDA < N c nc if NC > ncmax c c NOTE: the last column of JAO is used as work space!! c----------------------------------------------------------------------- integer i,j,k,ip real*8 zero logical copyval parameter (zero=0.0D0) c .. first executable statement .. copyval = (job.ne.0) if (lda .lt. n) then ierr = -1 return endif c .. use the last column of JAO as workspace c .. initialize the work space do i = 1, n jao(i,ncmax) = 0 enddo nc = 0 c .. go through ia and ja to find out number nonzero per row do k = 1, nnz i = ia(k) jao(i,ncmax) = jao(i,ncmax) + 1 enddo c .. maximum number of nonzero per row nc = 0 do i = 1, n if (nc.lt.jao(i,ncmax)) nc = jao(i,ncmax) jao(i,ncmax) = 0 enddo c .. if nc > ncmax retrun now if (nc.gt.ncmax) then ierr = nc return endif c .. go through ia and ja to copy the matrix to AO and JAO do k = 1, nnz i = ia(k) j = ja(k) jao(i,ncmax) = jao(i,ncmax) + 1 ip = jao(i,ncmax) if (ip.gt.nc) nc = ip if (copyval) ao(i,ip) = a(k) jao(i,ip) = j enddo c .. fill the unspecified elements of AO and JAO with zero diagonals do i = 1, n do j = ia(i+1)-ia(i)+1, nc jao(i,j)=i if(copyval) ao(i,j) = zero enddo enddo ierr = 0 c return end c-----end-of-cooell----------------------------------------------------- c----------------------------------------------------------------------- subroutine xcooell(n,nnz,a,ja,ia,ac,jac,nac,ner,ncmax,ierr) C----------------------------------------------------------------------- C coordinate format to ellpack format. C----------------------------------------------------------------------- C C DATE WRITTEN: June 4, 1989. C C PURPOSE C ------- C This subroutine takes a sparse matrix in coordinate format and C converts it into the Ellpack-Itpack storage. C C Example: C ------- C ( 11 0 13 0 0 0 ) C | 21 22 0 24 0 0 | C | 0 32 33 0 35 0 | C A = | 0 0 43 44 0 46 | C | 51 0 0 54 55 0 | C ( 61 62 0 0 65 66 ) C C Coordinate storage scheme: C C A = (11,22,33,44,55,66,13,21,24,32,35,43,46,51,54,61,62,65) C IA = (1, 2, 3, 4, 5, 6, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 6 ) C JA = ( 1, 2, 3, 4, 5, 6, 3, 1, 4, 2, 5, 3, 6, 1, 4, 1, 2, 5) C C Ellpack-Itpack storage scheme: C C ( 11 13 0 0 ) ( 1 3 * * ) C | 22 21 24 0 | | 2 1 4 * | C AC = | 33 32 35 0 | JAC = | 3 2 5 * | C | 44 43 46 0 | | 4 3 6 * | C | 55 51 54 0 | | 5 1 4 * | C ( 66 61 62 65 ) ( 6 1 2 5 ) C C Note: * means that you can store values from 1 to 6 (1 to n, where C n is the order of the matrix) in that position in the array. C C Contributed by: C --------------- C Ernest E. Rothman C Cornell Thoery Center/Cornell National Supercomputer Facility C e-mail address: BITNET: EER@CORNELLF.BITNET C INTERNET: eer@cornellf.tn.cornell.edu C C checked and modified 04/13/90 Y.Saad. C C REFERENCES C ---------- C Kincaid, D. R.; Oppe, T. C.; Respess, J. R.; Young, D. M. 1984. C ITPACKV 2C User's Guide, CNA-191. Center for Numerical Analysis, C University of Texas at Austin. C C "Engineering and Scientific Subroutine Library; Guide and C Reference; Release 3 (SC23-0184-3). Pp. 79-86. C C----------------------------------------------------------------------- C C INPUT PARAMETERS C ---------------- C N - Integer. The size of the square matrix. C C NNZ - Integer. Must be greater than or equal to the number of C nonzero elements in the sparse matrix. Dimension of A, IA C and JA. C C NCA - Integer. First dimension of output arrays ca and jac. C C A(NNZ) - Real array. (Double precision) C Stored entries of the sparse matrix A. C NNZ is the number of nonzeros. C C IA(NNZ) - Integer array. C Pointers to specify rows for the stored nonzero entries C in A. C C JA(NNZ) - Integer array. C Pointers to specify columns for the stored nonzero C entries in A. C C NER - Integer. Must be set greater than or equal to the maximum C number of nonzeros in any row of the sparse matrix. C C OUTPUT PARAMETERS C ----------------- C AC(NAC,*) - Real array. (Double precision) C Stored entries of the sparse matrix A in compressed C storage mode. C C JAC(NAC,*) - Integer array. C Contains the column numbers of the sparse matrix C elements stored in the corresponding positions in C array AC. C C NCMAX - Integer. Equals the maximum number of nonzeros in any C row of the sparse matrix. C C IERR - Error parameter is returned as zero on successful C execution of the subroutin