Skip to content

Commit 84e7123

Browse files
BoundaryValueProblemsBoundaryValueProblems
committed
Fixed diagonal matrix computations
1 parent 13fd242 commit 84e7123

File tree

2 files changed

+20
-16
lines changed

2 files changed

+20
-16
lines changed

‎src/HGLET.jl

Lines changed: 17 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -30,10 +30,10 @@ function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
3030
# Preliminaries
3131
W = G.W
3232
jmax = size(GP.rs,2)
33-
f = dmatrix[:,jmax,:]
3433

3534
# Fill in the appropriate entries of dmatrix
3635
dmatrix = dvec2dmatrix(dvec,GP,BS)
36+
f = dmatrix[:,jmax,:]
3737

3838
# Perform the signal synthesis from the given coefficients
3939
for j = jmax:-1:1
@@ -57,16 +57,16 @@ function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
5757
indrs = GP.ind[rs1:rs3-1]
5858
W_temp = W[indrs,indrs]
5959
D = Diagonal(vec(sum(W_temp, dims = 1)))
60-
D2 = D.^(-0.5)
60+
D2 =Diagonal(vec(sum(W_temp, dims = 1)).^(-0.5))
6161
normalizep = false # a flag for normalizing L for Lsym/Lrw
6262
if minimum(diag(D)) > 10^3*eps() # avoiding 1/small entries
6363
normalizep = true
6464
end
6565
# compute the GL eigenvectors via svd.
6666
if ( method == :Lrw || method == :Lsym ) && normalizep
67-
v,_,_ = svd(D2 * (D - W_temp) * D2)
67+
v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
6868
else
69-
v,_,_ = svd(D - W_temp)
69+
v,_,_ = svd(D - Matrix(W_temp))
7070
if method != :L
7171
@warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(method))
7272
end
@@ -82,6 +82,7 @@ function HGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
8282
standardized = true
8383
elseif v[row,col] < -10^3*eps()
8484
v[:,col] = - v[:,col]
85+
standardized = true
8586
else
8687
row += 1
8788
end
@@ -181,17 +182,17 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
181182
indrs = ind[rs1:rs3-1]
182183
W_temp = W[indrs,indrs]
183184
D = Diagonal(vec(sum(W_temp, dims = 1)))
184-
D2 = D.^(-0.5)
185+
D2 = Diagonal(vec(sum(W_temp, dims = 1)).^(-0.5))
185186
normalizep = false # a flag for normalizing L for Lsym and Lrw
186187
if minimum(diag(D)) > 10^3*eps() # avoiding 1/small entries
187188
normalizep = true
188189
end
189190

190191
# compute the GL eigenvectors via svd
191192
if ( method == :Lrw || method == :Lsym ) && normalizep
192-
v,_,_ = svd(D2 * (D - W_temp) * D2)
193+
v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
193194
else
194-
v,_,_ = svd(D - W_temp)
195+
v,_,_ = svd(D - Matrix(W_temp))
195196
if method != :L
196197
@warn("Due to the small diagonal entries of W, we revert back to the option :L, not :" * String(method))
197198
end
@@ -207,6 +208,7 @@ function HGLET_Analysis(G::GraphSig, GP::GraphPart, method::Symbol = :L)
207208
standardized = true
208209
elseif v[row,col] < -10^3*eps()
209210
v[:,col] = - v[:,col]
211+
standardized = true
210212
else
211213
row += 1
212214
end
@@ -274,7 +276,8 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
274276
indrs = ind[rs1:rs3-1]
275277

276278
# compute the eigenvectors of L ==> svd(L)
277-
v,_,_ = svd(Diagonal(vec(sum(W[indrs,indrs],dims = 1)))-W[indrs,indrs])
279+
v,_,_ = svd(Diagonal(vec(sum(W[indrs,indrs],dims = 1)))
280+
- Matrix(W[indrs,indrs]))
278281
v = v[:,end:-1:1]
279282

280283
# standardize the eigenvector signs
@@ -286,6 +289,7 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
286289
standardized = true
287290
elseif v[row,col] < -10^3*eps()
288291
v[:,col] = - v[:,col]
292+
standardized = true
289293
else
290294
row += 1
291295
end
@@ -325,17 +329,16 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
325329
### eigenvectors of L_rw ==> svd(L_sym)
326330
W_temp = W[indrs,indrs]
327331
D = Diagonal(vec(sum(W_temp,dims = 1)))
328-
D2 = D.^(-0.5)
329-
v,_,_ = svd(D2 * (D - W_temp) * D2)
332+
D2 =Diagonal(vec(sum(W_temp,dims = 1)).^(-0.5))
333+
v,_,_ = svd(D2 * (D - Matrix(W_temp)) * D2)
330334
v = v[:,end:-1:1]
331-
332335
else
333336
useLrw = false
334337

335338
### eigenvectors of L ==> svd(L)
336339
W_temp = W[indrs,indrs]
337340
D = Diagonal(vec(sum(W_temp,dims = 1)))
338-
v,_,_ = svd(D-W_temp)
341+
v,_,_ = svd(D - Matrix(W_temp))
339342
v = v[:,end:-1:1]
340343
end
341344

@@ -348,6 +351,7 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
348351
standardized = true
349352
elseif v[row,col] < -10^3*eps()
350353
v[:,col] = - v[:,col]
354+
standardized = true
351355
else
352356
row += 1
353357
end
@@ -359,7 +363,7 @@ function HGLET_Analysis_All(G::GraphSig, GP::GraphPart)
359363

360364
# obtain the expansion coeffcients for L_rw
361365
if useLrw
362-
dmatrixHrw[rs1:rs3-1,j,:] = v'*(Diagonal(vec(sum(W_temp,dims = 1))).^0.5)*G.f[indrs,:]
366+
dmatrixHrw[rs1:rs3-1,j,:] = v'*(Diagonal(vec(sum(W_temp,dims = 1)).^0.5))*G.f[indrs,:]
363367
else
364368
dmatrixHrw[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
365369
end

‎src/partition_fiedler.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ function partition_fiedler(W::SparseMatrixCSC{Float64,Int};
9191
end
9292
if N <= cutoff || eigs_flag != 0 # if W is small or eigs had a problem,
9393
# then use full svd.
94-
vtmp, val = svd(Diagonal(vec(sum(W, dims = 1))) - W) # v <-> U, val <-> S
94+
vtmp, val = svd(Matrix(Diagonal(vec(sum(W, dims = 1))) - W)) # v <-> U, val <-> S
9595
# Diagonal(vec(sum(W,1))) is Matrix{Float64} while W is SparseMatrix
9696
# But the results of the subtraction becomes Matrix{Float64}.
9797
# Also be careful here! val[end-2] == val[end-1] can happen.
@@ -133,8 +133,8 @@ function partition_fiedler(W::SparseMatrixCSC{Float64,Int};
133133
# then use full svd
134134
colsumW = vec(sum(W, dims = 1))
135135
D = Diagonal(colsumW)
136-
D2 =Diagonal(colsumW.^(-0.5))
137-
vtmp, val = svd(D2 * (D - W) * D2) # SVD of Lsym
136+
D2 = Diagonal(colsumW.^(-0.5))
137+
vtmp, val = svd(Matrix(D2 * (D - W) * D2)) # SVD of Lsym
138138
# MATLAB: [v,val,~] = svd(full(...
139139
# bsxfun(@times,bsxfun(@times,full(sum(W,2)).^(-0.5),diag(sum(W))-W),
140140
# full(sum(W,1)).^(-0.5)) ) );

0 commit comments

Comments
 (0)