Skip to content

Commit d5855d7

Browse files
mmuetzelmmuetzel
committed
mldivide: Correctly apply permutation for sparse matrix with permuted lower form.
* liboctave/array/CSparse.cc (SparseComplexMatrix::ltsolve), liboctave/array/dSparse.cc (SparseMatrix::ltsolve): Create and use inverse permutation vector for sparse matrix in permuted lower form. * test/mk-sparse-tst.sh: Add tests for mldivide with sparse matrices that have permuted lower form.
1 parent f648714 commit d5855d7

File tree

3 files changed

+149
-66
lines changed

3 files changed

+149
-66
lines changed

‎liboctave/array/CSparse.cc

Lines changed: 46 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -2588,14 +2588,18 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
25882588
{
25892589
retval.resize (nc, b_nc);
25902590
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2591+
25912592
octave_idx_type *perm = mattype.triangular_perm ();
2593+
OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
2594+
for (octave_idx_type i = 0; i < nr; i++)
2595+
inv_perm[perm[i]] = i;
25922596

25932597
for (octave_idx_type j = 0; j < b_nc; j++)
25942598
{
25952599
for (octave_idx_type i = 0; i < nm; i++)
25962600
work[i] = 0.;
25972601
for (octave_idx_type i = 0; i < nr; i++)
2598-
work[perm[i]] = b(i, j);
2602+
work[inv_perm[i]] = b(i, j);
25992603

26002604
for (octave_idx_type k = 0; k < nc; k++)
26012605
{
@@ -2604,10 +2608,11 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
26042608
octave_idx_type minr = nr;
26052609
octave_idx_type mini = 0;
26062610

2611+
// Check that permutation leads to lower triangular form
26072612
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2608-
if (perm[ridx (i)] < minr)
2613+
if (inv_perm[ridx (i)] < minr)
26092614
{
2610-
minr = perm[ridx (i)];
2615+
minr = inv_perm[ridx (i)];
26112616
mini = i;
26122617
}
26132618

@@ -2624,7 +2629,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
26242629
if (i == mini)
26252630
continue;
26262631

2627-
octave_idx_type iidx = perm[ridx (i)];
2632+
octave_idx_type iidx = inv_perm[ridx (i)];
26282633
work[iidx] = work[iidx] - tmp * data (i);
26292634
}
26302635
}
@@ -2653,9 +2658,9 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
26532658

26542659
for (octave_idx_type i = cidx (k);
26552660
i < cidx (k+1); i++)
2656-
if (perm[ridx (i)] < minr)
2661+
if (inv_perm[ridx (i)] < minr)
26572662
{
2658-
minr = perm[ridx (i)];
2663+
minr = inv_perm[ridx (i)];
26592664
mini = i;
26602665
}
26612666

@@ -2667,7 +2672,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
26672672
if (i == mini)
26682673
continue;
26692674

2670-
octave_idx_type iidx = perm[ridx (i)];
2675+
octave_idx_type iidx = inv_perm[ridx (i)];
26712676
work[iidx] = work[iidx] - tmp * data (i);
26722677
}
26732678
}
@@ -2841,14 +2846,18 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
28412846
if (typ == MatrixType::Permuted_Lower)
28422847
{
28432848
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2849+
28442850
octave_idx_type *perm = mattype.triangular_perm ();
2851+
OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
2852+
for (octave_idx_type i = 0; i < nr; i++)
2853+
inv_perm[perm[i]] = i;
28452854

28462855
for (octave_idx_type j = 0; j < b_nc; j++)
28472856
{
28482857
for (octave_idx_type i = 0; i < nm; i++)
28492858
work[i] = 0.;
28502859
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2851-
work[perm[b.ridx (i)]] = b.data (i);
2860+
work[inv_perm[b.ridx (i)]] = b.data (i);
28522861

28532862
for (octave_idx_type k = 0; k < nc; k++)
28542863
{
@@ -2857,10 +2866,11 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
28572866
octave_idx_type minr = nr;
28582867
octave_idx_type mini = 0;
28592868

2869+
// Check that permutation leads to lower triangular form
28602870
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2861-
if (perm[ridx (i)] < minr)
2871+
if (inv_perm[ridx (i)] < minr)
28622872
{
2863-
minr = perm[ridx (i)];
2873+
minr = inv_perm[ridx (i)];
28642874
mini = i;
28652875
}
28662876

@@ -2877,7 +2887,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
28772887
if (i == mini)
28782888
continue;
28792889

2880-
octave_idx_type iidx = perm[ridx (i)];
2890+
octave_idx_type iidx = inv_perm[ridx (i)];
28812891
work[iidx] = work[iidx] - tmp * data (i);
28822892
}
28832893
}
@@ -2928,9 +2938,9 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
29282938

29292939
for (octave_idx_type i = cidx (k);
29302940
i < cidx (k+1); i++)
2931-
if (perm[ridx (i)] < minr)
2941+
if (inv_perm[ridx (i)] < minr)
29322942
{
2933-
minr = perm[ridx (i)];
2943+
minr = inv_perm[ridx (i)];
29342944
mini = i;
29352945
}
29362946

@@ -2942,7 +2952,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
29422952
if (i == mini)
29432953
continue;
29442954

2945-
octave_idx_type iidx = perm[ridx (i)];
2955+
octave_idx_type iidx = inv_perm[ridx (i)];
29462956
work[iidx] = work[iidx] - tmp * data (i);
29472957
}
29482958
}
@@ -3134,14 +3144,18 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
31343144
{
31353145
retval.resize (nc, b_nc);
31363146
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3147+
31373148
octave_idx_type *perm = mattype.triangular_perm ();
3149+
OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
3150+
for (octave_idx_type i = 0; i < nr; i++)
3151+
inv_perm[perm[i]] = i;
31383152

31393153
for (octave_idx_type j = 0; j < b_nc; j++)
31403154
{
31413155
for (octave_idx_type i = 0; i < nm; i++)
31423156
work[i] = 0.;
31433157
for (octave_idx_type i = 0; i < nr; i++)
3144-
work[perm[i]] = b(i, j);
3158+
work[inv_perm[i]] = b(i, j);
31453159

31463160
for (octave_idx_type k = 0; k < nc; k++)
31473161
{
@@ -3151,9 +3165,9 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
31513165
octave_idx_type mini = 0;
31523166

31533167
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3154-
if (perm[ridx (i)] < minr)
3168+
if (inv_perm[ridx (i)] < minr)
31553169
{
3156-
minr = perm[ridx (i)];
3170+
minr = inv_perm[ridx (i)];
31573171
mini = i;
31583172
}
31593173

@@ -3170,7 +3184,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
31703184
if (i == mini)
31713185
continue;
31723186

3173-
octave_idx_type iidx = perm[ridx (i)];
3187+
octave_idx_type iidx = inv_perm[ridx (i)];
31743188
work[iidx] = work[iidx] - tmp * data (i);
31753189
}
31763190
}
@@ -3199,9 +3213,9 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
31993213

32003214
for (octave_idx_type i = cidx (k);
32013215
i < cidx (k+1); i++)
3202-
if (perm[ridx (i)] < minr)
3216+
if (inv_perm[ridx (i)] < minr)
32033217
{
3204-
minr = perm[ridx (i)];
3218+
minr = inv_perm[ridx (i)];
32053219
mini = i;
32063220
}
32073221

@@ -3213,7 +3227,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
32133227
if (i == mini)
32143228
continue;
32153229

3216-
octave_idx_type iidx = perm[ridx (i)];
3230+
octave_idx_type iidx = inv_perm[ridx (i)];
32173231
work[iidx] = work[iidx] - tmp * data (i);
32183232
}
32193233
}
@@ -3389,14 +3403,18 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
33893403
if (typ == MatrixType::Permuted_Lower)
33903404
{
33913405
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3406+
33923407
octave_idx_type *perm = mattype.triangular_perm ();
3408+
OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
3409+
for (octave_idx_type i = 0; i < nr; i++)
3410+
inv_perm[perm[i]] = i;
33933411

33943412
for (octave_idx_type j = 0; j < b_nc; j++)
33953413
{
33963414
for (octave_idx_type i = 0; i < nm; i++)
33973415
work[i] = 0.;
33983416
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3399-
work[perm[b.ridx (i)]] = b.data (i);
3417+
work[inv_perm[b.ridx (i)]] = b.data (i);
34003418

34013419
for (octave_idx_type k = 0; k < nc; k++)
34023420
{
@@ -3406,9 +3424,9 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
34063424
octave_idx_type mini = 0;
34073425

34083426
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3409-
if (perm[ridx (i)] < minr)
3427+
if (inv_perm[ridx (i)] < minr)
34103428
{
3411-
minr = perm[ridx (i)];
3429+
minr = inv_perm[ridx (i)];
34123430
mini = i;
34133431
}
34143432

@@ -3425,7 +3443,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
34253443
if (i == mini)
34263444
continue;
34273445

3428-
octave_idx_type iidx = perm[ridx (i)];
3446+
octave_idx_type iidx = inv_perm[ridx (i)];
34293447
work[iidx] = work[iidx] - tmp * data (i);
34303448
}
34313449
}
@@ -3476,9 +3494,9 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
34763494

34773495
for (octave_idx_type i = cidx (k);
34783496
i < cidx (k+1); i++)
3479-
if (perm[ridx (i)] < minr)
3497+
if (inv_perm[ridx (i)] < minr)
34803498
{
3481-
minr = perm[ridx (i)];
3499+
minr = inv_perm[ridx (i)];
34823500
mini = i;
34833501
}
34843502

@@ -3490,7 +3508,7 @@ SparseComplexMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
34903508
if (i == mini)
34913509
continue;
34923510

3493-
octave_idx_type iidx = perm[ridx (i)];
3511+
octave_idx_type iidx = inv_perm[ridx (i)];
34943512
work[iidx] = work[iidx] - tmp * data (i);
34953513
}
34963514
}

0 commit comments

Comments
 (0)