@@ -104,11 +104,26 @@ struct BidiagonalIndex <: MatrixIndex
104
104
end
105
105
106
106
struct TridiagonalIndex <: MatrixIndex
107
+ count:: Int # count==nsize+nsize-1+nsize-1
108
+ nsize:: Int
109
+ isrow:: Bool
110
+ end
111
+
112
+ struct BandedMatrixIndex <: MatrixIndex
107
113
count:: Int
108
114
nsize:: Int
115
+ bandinds:: Array{Int}
116
+ bandsizes:: Array{Int}
109
117
isrow:: Bool
110
118
end
111
119
120
+ function BandedMatrixIndex (nsize,lowerbandwidth,upperbandwidth,isrow)
121
+ upperbandwidth> - lowerbandwidth || throw (ErrorException (" Invalid Bandwidths" ))
122
+ bandinds= upperbandwidth: - 1 : - lowerbandwidth
123
+ bandsizes= [nsize- abs (band) for band in bandinds]
124
+ BandedMatrixIndex (sum (bandsizes),nsize,bandinds,bandsizes,isrow)
125
+ end
126
+
112
127
Base. firstindex (ind:: MatrixIndex )= 1
113
128
Base. lastindex (ind:: MatrixIndex )= ind. count
114
129
Base. length (ind:: MatrixIndex )= ind. count
@@ -135,6 +150,23 @@ function Base.getindex(ind::TridiagonalIndex,i::Int)
135
150
end
136
151
end
137
152
153
+ function Base. getindex (ind:: BandedMatrixIndex ,i:: Int )
154
+ 1 <= i <= ind. count || throw (BoundsError (ind, i))
155
+ _i= i
156
+ p= 1
157
+ while _i- ind. bandsizes[p]> 0
158
+ _i-= ind. bandsizes[p]
159
+ p+= 1
160
+ end
161
+ bandind= ind. bandinds[p]
162
+ startfromone= ind. isrow & (bandind> 0 )
163
+ if startfromone
164
+ return _i
165
+ else
166
+ return _i+ abs (bandind)
167
+ end
168
+ end
169
+
138
170
function findstructralnz (x:: Bidiagonal )
139
171
n= size (x,1 )
140
172
isup= x. uplo== ' U' ? true : false
@@ -225,11 +257,20 @@ function __init__()
225
257
end
226
258
227
259
@require BandedMatrices= " aae01518-5342-5314-be14-df237901396f" begin
260
+ function findstructralnz (x:: BandedMatrices.BandedMatrix )
261
+ l,u= BandedMatrices. bandwidths (x)
262
+ nsize= size (x,1 )
263
+ rowind= BandedMatrixIndex (nsize,l,u,true )
264
+ colind= BandedMatrixIndex (nsize,l,u,false )
265
+ (rowind,colind)
266
+ end
267
+
268
+ has_sparsestruct (:: Type{<:BandedMatrices.BandedMatrix} ) = true
228
269
is_structured (:: Type{<:BandedMatrices.BandedMatrix} ) = true
229
270
fast_matrix_colors (:: Type{<:BandedMatrices.BandedMatrix} ) = true
230
271
231
272
function matrix_colors (A:: BandedMatrices.BandedMatrix )
232
- u,l = bandwidths (A)
273
+ l,u = BandedMatrices . bandwidths (A)
233
274
width= u+ l+ 1
234
275
_cycle (1 : width,size (A,2 ))
235
276
end
@@ -243,10 +284,10 @@ function __init__()
243
284
fast_matrix_colors (:: Type{<:BlockBandedMatrices.BandedBlockBandedMatrix} ) = true
244
285
245
286
function matrix_colors (A:: BlockBandedMatrices.BlockBandedMatrix )
246
- l,u= blockbandwidths (A)
287
+ l,u= BlockBandedMatrices . blockbandwidths (A)
247
288
blockwidth= l+ u+ 1
248
- nblock= nblocks (A,2 )
249
- cols= [blocksize (A,(1 ,i))[2 ] for i in 1 : nblock]
289
+ nblock= BlockBandedMatrices . nblocks (A,2 )
290
+ cols= [BlockBandedMatrices . blocksize (A,(1 ,i))[2 ] for i in 1 : nblock]
250
291
blockcolors= _cycle (1 : blockwidth,nblock)
251
292
# the reserved number of colors of a block is the maximum length of columns of blocks with the same block color
252
293
ncolors= [maximum (cols[i: blockwidth: nblock]) for i in 1 : blockwidth]
@@ -257,12 +298,12 @@ function __init__()
257
298
end
258
299
259
300
function matrix_colors (A:: BlockBandedMatrices.BandedBlockBandedMatrix )
260
- l,u= blockbandwidths (A)
261
- lambda,mu= subblockbandwidths (A)
301
+ l,u= BlockBandedMatrices . blockbandwidths (A)
302
+ lambda,mu= BlockBandedMatrices . subblockbandwidths (A)
262
303
blockwidth= l+ u+ 1
263
304
subblockwidth= lambda+ mu+ 1
264
- nblock= nblocks (A,2 )
265
- cols= [blocksize (A,(1 ,i))[2 ] for i in 1 : nblock]
305
+ nblock= BlockBandedMatrices . nblocks (A,2 )
306
+ cols= [BlockBandedMatrices . blocksize (A,(1 ,i))[2 ] for i in 1 : nblock]
266
307
blockcolors= _cycle (1 : blockwidth,nblock)
267
308
# the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color
268
309
ncolors= [min (subblockwidth,maximum (cols[i: blockwidth: nblock])) for i in 1 : min (blockwidth,nblock)]
0 commit comments