@@ -805,6 +805,53 @@ function __init__()
805805 end
806806
807807 @require BandedMatrices = " aae01518-5342-5314-be14-df237901396f" begin
808+ struct BandedMatrixIndex <: MatrixIndex
809+ count:: Int
810+ rowsize:: Int
811+ colsize:: Int
812+ bandinds:: Array{Int,1}
813+ bandsizes:: Array{Int,1}
814+ isrow:: Bool
815+ end
816+
817+ Base. firstindex (i:: BandedMatrixIndex ) = 1
818+ Base. lastindex (i:: BandedMatrixIndex ) = i. count
819+ Base. length (i:: BandedMatrixIndex ) = lastindex (i)
820+ @propagate_inbounds function Base. getindex (ind:: BandedMatrixIndex , i:: Int )
821+ @boundscheck 1 <= i <= ind. count || throw (BoundsError (ind, i))
822+ _i = i
823+ p = 1
824+ while _i - ind. bandsizes[p] > 0
825+ _i -= ind. bandsizes[p]
826+ p += 1
827+ end
828+ bandind = ind. bandinds[p]
829+ startfromone = ! xor (ind. isrow, (bandind > 0 ))
830+ if startfromone
831+ return _i
832+ else
833+ return _i + abs (bandind)
834+ end
835+ end
836+
837+ function _bandsize (bandind, rowsize, colsize)
838+ - (rowsize - 1 ) <= bandind <= colsize - 1 || throw (ErrorException (" Invalid Bandind" ))
839+ if (bandind * (colsize - rowsize) > 0 ) & (abs (bandind) <= abs (colsize - rowsize))
840+ return min (rowsize, colsize)
841+ elseif bandind * (colsize - rowsize) <= 0
842+ return min (rowsize, colsize) - abs (bandind)
843+ else
844+ return min (rowsize, colsize) - abs (bandind) + abs (colsize - rowsize)
845+ end
846+ end
847+
848+ function BandedMatrixIndex (rowsize, colsize, lowerbandwidth, upperbandwidth, isrow)
849+ upperbandwidth > - lowerbandwidth || throw (ErrorException (" Invalid Bandwidths" ))
850+ bandinds = upperbandwidth: - 1 : - lowerbandwidth
851+ bandsizes = [_bandsize (band, rowsize, colsize) for band in bandinds]
852+ BandedMatrixIndex (sum (bandsizes), rowsize, colsize, bandinds, bandsizes, isrow)
853+ end
854+
808855 function findstructralnz (x:: BandedMatrices.BandedMatrix )
809856 l, u = BandedMatrices. bandwidths (x)
810857 rowsize, colsize = Base. size (x)
@@ -827,6 +874,82 @@ function __init__()
827874
828875 @require BlockBandedMatrices = " ffab5731-97b5-5995-9138-79e8c1846df0" begin
829876 @require BlockArrays = " 8e7c35d0-a365-5155-bbbb-fb81a777f24e" begin
877+ struct BlockBandedMatrixIndex <: MatrixIndex
878+ count:: Int
879+ refinds:: Array{Int,1}
880+ refcoords:: Array{Int,1} # storing col or row inds at ref points
881+ isrow:: Bool
882+ end
883+ Base. firstindex (i:: BlockBandedMatrixIndex ) = 1
884+ Base. lastindex (i:: BlockBandedMatrixIndex ) = i. count
885+ Base. length (i:: BlockBandedMatrixIndex ) = lastindex (i)
886+ function BlockBandedMatrixIndex (nrowblock, ncolblock, rowsizes, colsizes, l, u)
887+ blockrowind = BandedMatrixIndex (nrowblock, ncolblock, l, u, true )
888+ blockcolind = BandedMatrixIndex (nrowblock, ncolblock, l, u, false )
889+ sortedinds = sort (
890+ [(blockrowind[i], blockcolind[i]) for i = 1 : length (blockrowind)],
891+ by = x -> x[1 ],
892+ )
893+ sort! (sortedinds, by = x -> x[2 ], alg = InsertionSort)# stable sort keeps the second index in order
894+ refinds = Array {Int,1} ()
895+ refrowcoords = Array {Int,1} ()
896+ refcolcoords = Array {Int,1} ()
897+ rowheights = pushfirst! (copy (rowsizes), 1 )
898+ cumsum! (rowheights, rowheights)
899+ blockheight = 0
900+ blockrow = 1
901+ blockcol = 1
902+ currenti = 1
903+ lastrowind = sortedinds[1 ][1 ] - 1
904+ lastcolind = sortedinds[1 ][2 ]
905+ for ind in sortedinds
906+ rowind, colind = ind
907+ if colind == lastcolind
908+ if rowind > lastrowind
909+ blockheight += rowsizes[rowind]
910+ end
911+ else
912+ for j = blockcol: blockcol+ colsizes[lastcolind]- 1
913+ push! (refinds, currenti)
914+ push! (refrowcoords, blockrow)
915+ push! (refcolcoords, j)
916+ currenti += blockheight
917+ end
918+ blockcol += colsizes[lastcolind]
919+ blockrow = rowheights[rowind]
920+ blockheight = rowsizes[rowind]
921+ end
922+ lastcolind = colind
923+ lastrowind = rowind
924+ end
925+ for j = blockcol: blockcol+ colsizes[lastcolind]- 1
926+ push! (refinds, currenti)
927+ push! (refrowcoords, blockrow)
928+ push! (refcolcoords, j)
929+ currenti += blockheight
930+ end
931+ push! (refinds, currenti)# guard
932+ push! (refrowcoords, - 1 )
933+ push! (refcolcoords, - 1 )
934+ rowindobj = BlockBandedMatrixIndex (currenti - 1 , refinds, refrowcoords, true )
935+ colindobj = BlockBandedMatrixIndex (currenti - 1 , refinds, refcolcoords, false )
936+ rowindobj, colindobj
937+ end
938+ @propagate_inbounds function Base. getindex (ind:: BlockBandedMatrixIndex , i:: Int )
939+ @boundscheck 1 <= i <= ind. count || throw (BoundsError (ind, i))
940+ p = 1
941+ while i - ind. refinds[p] >= 0
942+ p += 1
943+ end
944+ p -= 1
945+ _i = i - ind. refinds[p]
946+ if ind. isrow
947+ return ind. refcoords[p] + _i
948+ else
949+ return ind. refcoords[p]
950+ end
951+ end
952+
830953 function findstructralnz (x:: BlockBandedMatrices.BlockBandedMatrix )
831954 l, u = BlockBandedMatrices. blockbandwidths (x)
832955 nrowblock = BlockBandedMatrices. blocksize (x, 1 )
@@ -842,6 +965,87 @@ function __init__()
842965 u,
843966 )
844967 end
968+ struct BandedBlockBandedMatrixIndex <: MatrixIndex
969+ count:: Int
970+ refinds:: Array{Int,1}
971+ refcoords:: Array{Int,1} # storing col or row inds at ref points
972+ reflocalinds:: Array{BandedMatrixIndex,1}
973+ isrow:: Bool
974+ end
975+ Base. firstindex (i:: BandedBlockBandedMatrixIndex ) = 1
976+ Base. lastindex (i:: BandedBlockBandedMatrixIndex ) = i. count
977+ Base. length (i:: BandedBlockBandedMatrixIndex ) = lastindex (i)
978+ @propagate_inbounds function Base. getindex (ind:: BandedBlockBandedMatrixIndex , i:: Int )
979+ @boundscheck 1 <= i <= ind. count || throw (BoundsError (ind, i))
980+ p = 1
981+ while i - ind. refinds[p] >= 0
982+ p += 1
983+ end
984+ p -= 1
985+ _i = i - ind. refinds[p] + 1
986+ ind. reflocalinds[p][_i] + ind. refcoords[p] - 1
987+ end
988+
989+
990+ function BandedBlockBandedMatrixIndex (
991+ nrowblock,
992+ ncolblock,
993+ rowsizes,
994+ colsizes,
995+ l,
996+ u,
997+ lambda,
998+ mu,
999+ )
1000+ blockrowind = BandedMatrixIndex (nrowblock, ncolblock, l, u, true )
1001+ blockcolind = BandedMatrixIndex (nrowblock, ncolblock, l, u, false )
1002+ sortedinds = sort (
1003+ [(blockrowind[i], blockcolind[i]) for i = 1 : length (blockrowind)],
1004+ by = x -> x[1 ],
1005+ )
1006+ sort! (sortedinds, by = x -> x[2 ], alg = InsertionSort)# stable sort keeps the second index in order
1007+ rowheights = pushfirst! (copy (rowsizes), 1 )
1008+ cumsum! (rowheights, rowheights)
1009+ colwidths = pushfirst! (copy (colsizes), 1 )
1010+ cumsum! (colwidths, colwidths)
1011+ currenti = 1
1012+ refinds = Array {Int,1} ()
1013+ refrowcoords = Array {Int,1} ()
1014+ refcolcoords = Array {Int,1} ()
1015+ reflocalrowinds = Array {BandedMatrixIndex,1} ()
1016+ reflocalcolinds = Array {BandedMatrixIndex,1} ()
1017+ for ind in sortedinds
1018+ rowind, colind = ind
1019+ localrowind =
1020+ BandedMatrixIndex (rowsizes[rowind], colsizes[colind], lambda, mu, true )
1021+ localcolind =
1022+ BandedMatrixIndex (rowsizes[rowind], colsizes[colind], lambda, mu, false )
1023+ push! (refinds, currenti)
1024+ push! (refrowcoords, rowheights[rowind])
1025+ push! (refcolcoords, colwidths[colind])
1026+ push! (reflocalrowinds, localrowind)
1027+ push! (reflocalcolinds, localcolind)
1028+ currenti += localrowind. count
1029+ end
1030+ push! (refinds, currenti)
1031+ push! (refrowcoords, - 1 )
1032+ push! (refcolcoords, - 1 )
1033+ rowindobj = BandedBlockBandedMatrixIndex (
1034+ currenti - 1 ,
1035+ refinds,
1036+ refrowcoords,
1037+ reflocalrowinds,
1038+ true ,
1039+ )
1040+ colindobj = BandedBlockBandedMatrixIndex (
1041+ currenti - 1 ,
1042+ refinds,
1043+ refcolcoords,
1044+ reflocalcolinds,
1045+ false ,
1046+ )
1047+ rowindobj, colindobj
1048+ end
8451049
8461050 function findstructralnz (x:: BlockBandedMatrices.BandedBlockBandedMatrix )
8471051 l, u = BlockBandedMatrices. blockbandwidths (x)
0 commit comments