1 /*
2 * Copyright © 2004-2011 Ondra Kamenik
3 * Copyright © 2019 Dynare Team
4 *
5 * This file is part of Dynare.
6 *
7 * Dynare is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Dynare is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include "QuasiTriangular.hh"
22 #include "SylvException.hh"
23 #include "SchurDecomp.hh"
24 #include "int_power.hh"
25
26 #include <dynblas.h>
27
28 #include <cmath>
29 #include <iostream>
30 #include <sstream>
31
32 double
getDeterminant() const33 DiagonalBlock::getDeterminant() const
34 {
35 return (*alpha)*(*alpha) + getSBeta();
36 }
37
38 double
getSBeta() const39 DiagonalBlock::getSBeta() const
40 {
41 return -(*beta1)*(*beta2);
42 }
43
44 double
getSize() const45 DiagonalBlock::getSize() const
46 {
47 if (real)
48 return std::abs(*alpha);
49 else
50 return std::sqrt(getDeterminant());
51 }
52
53 /* This function makes Diagonal inconsistent, it should only be used
54 on temorary matrices, which will not be used any more, e.g. in
55 QuasiTriangular::solve() (we need fast performance) */
56 void
setReal()57 DiagonalBlock::setReal()
58 {
59 *beta1 = 0;
60 *beta2 = 0;
61 real = true;
62 }
63
64 void
checkBlock(const double * d,int d_size)65 DiagonalBlock::checkBlock(const double *d, int d_size)
66 {
67 const double *a1 = d + jbar*d_size+jbar;
68 const double *b1 = a1 + d_size;
69 const double *b2 = a1 + 1;
70 const double *a2 = b1 + 1;
71 if (a1 != alpha.a1)
72 throw SYLV_MES_EXCEPTION("Bad alpha1.");
73 if (!real && b1 != beta1)
74 throw SYLV_MES_EXCEPTION("Bad beta1.");
75 if (!real && b2 != beta2)
76 throw SYLV_MES_EXCEPTION("Bad beta2.");
77 if (!real && a2 != alpha.a2)
78 throw SYLV_MES_EXCEPTION("Bad alpha2.");
79 }
80
Diagonal(double * data,int d_size)81 Diagonal::Diagonal(double *data, int d_size)
82 {
83 int nc = getNumComplex(data, d_size); // return nc ≤ d_size/2
84 num_all = d_size - nc;
85 num_real = d_size - 2*nc;
86
87 int jbar = 0;
88 int j = 0;
89 while (j < num_all)
90 {
91 int id = jbar*d_size + jbar; // index of diagonal block in data
92 int ill = id + 1; // index of element below the diagonal
93 int iur = id + d_size; // index of element right to diagonal
94 int idd = id + d_size + 1; // index of element next on diagonal
95 if ((jbar < d_size-1) && !isZero(data[ill]))
96 {
97 // it is not last column and we have nonzero below diagonal
98 blocks.emplace_back(jbar, false, &data[id], &data[idd],
99 &data[iur], &data[ill]);
100 jbar++;
101 }
102 else
103 // it is last column or we have zero below diagonal
104 blocks.emplace_back(jbar, true, &data[id], &data[id], nullptr, nullptr);
105 jbar++;
106 j++;
107 }
108 }
109
Diagonal(double * data,const Diagonal & d)110 Diagonal::Diagonal(double *data, const Diagonal &d)
111 {
112 num_all = d.num_all;
113 num_real = d.num_real;
114 int d_size = d.getSize();
115 for (const auto &block : d)
116 {
117 double *beta1 = nullptr;
118 double *beta2 = nullptr;
119 int id = block.getIndex()*(d_size+1);
120 int idd = id;
121 if (!block.isReal())
122 {
123 beta1 = &data[id+d_size];
124 beta2 = &data[id+1];
125 idd = id + d_size + 1;
126 }
127 blocks.emplace_back(block.getIndex(), block.isReal(),
128 &data[id], &data[idd], beta1, beta2);
129 }
130 }
131
132 int
getNumComplex(const double * data,int d_size)133 Diagonal::getNumComplex(const double *data, int d_size)
134 {
135 int num_complex = 0;
136 int in = 1;
137 for (int i = 0; i < d_size-1; i++, in = in + d_size + 1)
138 if (!isZero(data[in]))
139 {
140 num_complex++;
141 if (in < d_size - 2 && !isZero(data[in + d_size +1]))
142 throw SYLV_MES_EXCEPTION("Matrix is not quasi-triangular");
143 }
144 return num_complex;
145 }
146
147 void
changeBase(double * p)148 Diagonal::changeBase(double *p)
149 {
150 int d_size = getSize();
151 for (auto &it : *this)
152 {
153 const DiagonalBlock &b = it;
154 int jbar = b.getIndex();
155 int base = d_size*jbar + jbar;
156 if (b.isReal())
157 {
158 DiagonalBlock bnew(jbar, true, &p[base], &p[base],
159 nullptr, nullptr);
160 it = bnew;
161 }
162 else
163 {
164 DiagonalBlock bnew(jbar, false, &p[base], &p[base+d_size+1],
165 &p[base+d_size], &p[base+1]);
166 it = bnew;
167 }
168 }
169 }
170
171 void
getEigenValues(Vector & eig) const172 Diagonal::getEigenValues(Vector &eig) const
173 {
174 int d_size = getSize();
175 if (eig.length() != 2*d_size)
176 {
177 std::ostringstream mes;
178 mes << "Wrong length of vector for eigenvalues len=" << eig.length()
179 << ", should be=" << 2*d_size << '.' << std::endl;
180 throw SYLV_MES_EXCEPTION(mes.str());
181 }
182 for (const auto &b : *this)
183 {
184 int ind = b.getIndex();
185 eig[2*ind] = *(b.getAlpha());
186 if (b.isReal())
187 eig[2*ind+1] = 0.0;
188 else
189 {
190 double beta = std::sqrt(b.getSBeta());
191 eig[2*ind+1] = beta;
192 eig[2*ind+2] = eig[2*ind];
193 eig[2*ind+3] = -beta;
194 }
195 }
196 }
197
198 /* Swaps logically blocks ‘it’, and ‘++it’. remember to move also
199 addresses, alpha, beta1, beta2. This is a dirty (but most
200 effective) way how to do it. */
201 void
swapLogically(diag_iter it)202 Diagonal::swapLogically(diag_iter it)
203 {
204 diag_iter itp = it;
205 ++itp;
206
207 if (it->isReal() && !itp->isReal())
208 {
209 // first is real, second is complex
210 double *d1 = it->alpha.a1;
211 double *d2 = itp->alpha.a1;
212 double *d3 = itp->alpha.a2;
213 // swap
214 DiagonalBlock new_it(it->jbar, d1, d2);
215 *it = new_it;
216 DiagonalBlock new_itp(itp->jbar+1, d3);
217 *itp = new_itp;
218 }
219 else if (!it->isReal() && itp->isReal())
220 {
221 // first is complex, second is real
222 double *d1 = it->alpha.a1;
223 double *d2 = it->alpha.a2;
224 double *d3 = itp->alpha.a1;
225 // swap
226 DiagonalBlock new_it(it->jbar, d1);
227 *it = new_it;
228 DiagonalBlock new_itp(itp->jbar-1, d2, d3);
229 *itp = new_itp;
230 }
231 }
232
233 void
checkConsistency(diag_iter it)234 Diagonal::checkConsistency(diag_iter it)
235 {
236 if (!it->isReal() && isZero(it->getBeta2()))
237 {
238 it->getBeta2() = 0.0; // put exact zero
239 int jbar = it->getIndex();
240 double *d2 = it->alpha.a2;
241 it->alpha.a2 = it->alpha.a1;
242 it->real = true;
243 it->beta1 = nullptr;
244 it->beta2 = nullptr;
245 blocks.emplace(++it, jbar+1, d2);
246 num_real += 2;
247 num_all++;
248 }
249 }
250
251 double
getAverageSize(diag_iter start,diag_iter end)252 Diagonal::getAverageSize(diag_iter start, diag_iter end)
253 {
254 double res = 0;
255 int num = 0;
256 for (diag_iter run = start; run != end; ++run)
257 {
258 num++;
259 res += run->getSize();
260 }
261 if (num > 0)
262 res = res/num;
263 return res;
264 }
265
266 Diagonal::diag_iter
findClosestBlock(diag_iter start,diag_iter end,double a)267 Diagonal::findClosestBlock(diag_iter start, diag_iter end, double a)
268 {
269 diag_iter closest = start;
270 double minim = 1.0e100;
271 for (diag_iter run = start; run != end; ++run)
272 {
273 double dist = std::abs(a - run->getSize());
274 if (dist < minim)
275 {
276 minim = dist;
277 closest = run;
278 }
279 }
280 return closest;
281 }
282
283 Diagonal::diag_iter
findNextLargerBlock(diag_iter start,diag_iter end,double a)284 Diagonal::findNextLargerBlock(diag_iter start, diag_iter end, double a)
285 {
286 diag_iter closest = start;
287 double minim = 1.0e100;
288 for (diag_iter run = start; run != end; ++run)
289 {
290 double dist = run->getSize() - a;
291 if ((0 <= dist) && (dist < minim))
292 {
293 minim = dist;
294 closest = run;
295 }
296 }
297 return closest;
298 }
299
300 void
print() const301 Diagonal::print() const
302 {
303 auto ff = std::cout.flags();
304 std::cout << "Num real: " << getNumReal() << ", num complex: " << getNumComplex() << std::endl
305 << std::fixed;
306 for (const auto &it : *this)
307 if (it.isReal())
308 std::cout << "real: jbar=" << it.getIndex() << ", alpha=" << *(it.getAlpha()) << std::endl;
309 else
310 std::cout << "complex: jbar=" << it.getIndex()
311 << ", alpha=" << *(it.getAlpha())
312 << ", beta1=" << it.getBeta1()
313 << ", beta2=" << it.getBeta2() << std::endl;
314 std::cout.flags(ff);
315 }
316
317 bool
isZero(double p)318 Diagonal::isZero(double p)
319 {
320 return (std::abs(p) < EPS);
321 }
322
323 QuasiTriangular::const_col_iter
col_begin(const DiagonalBlock & b) const324 QuasiTriangular::col_begin(const DiagonalBlock &b) const
325 {
326 int jbar = b.getIndex();
327 int d_size = diagonal.getSize();
328 return const_col_iter(&getData()[jbar*d_size], d_size, b.isReal(), 0);
329 }
330
331 QuasiTriangular::col_iter
col_begin(const DiagonalBlock & b)332 QuasiTriangular::col_begin(const DiagonalBlock &b)
333 {
334 int jbar = b.getIndex();
335 int d_size = diagonal.getSize();
336 return col_iter(&getData()[jbar*d_size], d_size, b.isReal(), 0);
337 }
338
339 QuasiTriangular::const_row_iter
row_begin(const DiagonalBlock & b) const340 QuasiTriangular::row_begin(const DiagonalBlock &b) const
341 {
342 int jbar = b.getIndex();
343 int d_size = diagonal.getSize();
344 int off = jbar*d_size+jbar+d_size;
345 int col = jbar+1;
346 if (!b.isReal())
347 {
348 off = off + d_size;
349 col++;
350 }
351 return const_row_iter(&getData()[off], d_size, b.isReal(), col);
352 }
353
354 QuasiTriangular::row_iter
row_begin(const DiagonalBlock & b)355 QuasiTriangular::row_begin(const DiagonalBlock &b)
356 {
357 int jbar = b.getIndex();
358 int d_size = diagonal.getSize();
359 int off = jbar*d_size+jbar+d_size;
360 int col = jbar+1;
361 if (!b.isReal())
362 {
363 off = off + d_size;
364 col++;
365 }
366 return row_iter(&getData()[off], d_size, b.isReal(), col);
367 }
368
369 QuasiTriangular::const_col_iter
col_end(const DiagonalBlock & b) const370 QuasiTriangular::col_end(const DiagonalBlock &b) const
371 {
372 int jbar = b.getIndex();
373 int d_size = diagonal.getSize();
374 return const_col_iter(getData().base()+jbar*d_size+jbar, d_size, b.isReal(),
375 jbar);
376 }
377
378 QuasiTriangular::col_iter
col_end(const DiagonalBlock & b)379 QuasiTriangular::col_end(const DiagonalBlock &b)
380 {
381 int jbar = b.getIndex();
382 int d_size = diagonal.getSize();
383 return col_iter(&getData()[jbar*d_size+jbar], d_size, b.isReal(), jbar);
384 }
385
386 QuasiTriangular::const_row_iter
row_end(const DiagonalBlock & b) const387 QuasiTriangular::row_end(const DiagonalBlock &b) const
388 {
389 int jbar = b.getIndex();
390 int d_size = diagonal.getSize();
391 return const_row_iter(&getData()[d_size*d_size+jbar], d_size, b.isReal(),
392 d_size);
393 }
394
395 QuasiTriangular::row_iter
row_end(const DiagonalBlock & b)396 QuasiTriangular::row_end(const DiagonalBlock &b)
397 {
398 int jbar = b.getIndex();
399 int d_size = diagonal.getSize();
400 return row_iter(&getData()[d_size*d_size+jbar], d_size, b.isReal(), d_size);
401 }
402
QuasiTriangular(double r,const QuasiTriangular & t)403 QuasiTriangular::QuasiTriangular(double r, const QuasiTriangular &t)
404 : SqSylvMatrix(t.nrows()), diagonal(getData().base(), t.diagonal)
405 {
406 setMatrix(r, t);
407 }
408
QuasiTriangular(double r,const QuasiTriangular & t,double r2,const QuasiTriangular & t2)409 QuasiTriangular::QuasiTriangular(double r, const QuasiTriangular &t,
410 double r2, const QuasiTriangular &t2)
411 : SqSylvMatrix(t.nrows()), diagonal(getData().base(), t.diagonal)
412 {
413 setMatrix(r, t);
414 addMatrix(r2, t2);
415 }
416
QuasiTriangular(const QuasiTriangular & t)417 QuasiTriangular::QuasiTriangular(const QuasiTriangular &t)
418 : SqSylvMatrix(t), diagonal(getData().base(), t.diagonal)
419 {
420 }
421
QuasiTriangular(const ConstVector & d,int d_size)422 QuasiTriangular::QuasiTriangular(const ConstVector &d, int d_size)
423 : SqSylvMatrix(Vector{d}, d_size), diagonal(getData().base(), d_size)
424 {
425 }
426
QuasiTriangular(const std::string & dummy,const QuasiTriangular & t)427 QuasiTriangular::QuasiTriangular(const std::string &dummy, const QuasiTriangular &t)
428 : SqSylvMatrix(t.nrows()), diagonal(getData().base(), t.diagonal)
429 {
430 Vector aux(t.getData());
431 blas_int d_size = diagonal.getSize();
432 double alpha = 1.0;
433 double beta = 0.0;
434 blas_int lda = t.getLD(), ldb = t.getLD(), ldc = ld;
435 dgemm("N", "N", &d_size, &d_size, &d_size, &alpha, aux.base(),
436 &lda, t.getData().base(), &ldb, &beta, getData().base(), &ldc);
437 }
438
QuasiTriangular(const SchurDecomp & decomp)439 QuasiTriangular::QuasiTriangular(const SchurDecomp &decomp)
440 : SqSylvMatrix(decomp.getT()),
441 diagonal(getData().base(), decomp.getDim())
442 {
443 }
444
445 // This pads matrix with intial columns with zeros
QuasiTriangular(const SchurDecompZero & decomp)446 QuasiTriangular::QuasiTriangular(const SchurDecompZero &decomp)
447 : SqSylvMatrix(decomp.getDim())
448 {
449 // nullify first decomp.getZeroCols() columns
450 int zeros = decomp.getZeroCols()*decomp.getDim();
451 Vector zv(getData(), 0, zeros);
452 zv.zeros();
453 // fill right upper part with decomp.getRU()
454 for (int i = 0; i < decomp.getRU().nrows(); i++)
455 for (int j = 0; j < decomp.getRU().ncols(); j++)
456 getData()[(j+decomp.getZeroCols())*decomp.getDim()+i] = decomp.getRU().get(i, j);
457
458 // fill right lower part with decomp.getT()
459 for (int i = 0; i < decomp.getT().nrows(); i++)
460 for (int j = 0; j < decomp.getT().ncols(); j++)
461 getData()[(j+decomp.getZeroCols())*decomp.getDim()+decomp.getZeroCols()+i]
462 = decomp.getT().get(i, j);
463
464 // construct diagonal
465 diagonal = Diagonal{getData().base(), decomp.getDim()};
466 }
467
468 void
setMatrix(double r,const QuasiTriangular & t)469 QuasiTriangular::setMatrix(double r, const QuasiTriangular &t)
470 {
471 getData().zeros();
472 getData().add(r, t.getData());
473 }
474
475 void
addMatrix(double r,const QuasiTriangular & t)476 QuasiTriangular::addMatrix(double r, const QuasiTriangular &t)
477 {
478 getData().add(r, t.getData());
479 }
480
481 void
addUnit()482 QuasiTriangular::addUnit()
483 {
484 for (diag_iter di = diag_begin(); di != diag_end(); ++di)
485 di->getAlpha() = *(di->getAlpha()) + 1.0;
486 }
487
488 void
solve(Vector & x,const ConstVector & b,double & eig_min)489 QuasiTriangular::solve(Vector &x, const ConstVector &b, double &eig_min)
490 {
491 x = b;
492 solvePre(x, eig_min);
493 }
494
495 void
solveTrans(Vector & x,const ConstVector & b,double & eig_min)496 QuasiTriangular::solveTrans(Vector &x, const ConstVector &b, double &eig_min)
497 {
498 x = b;
499 solvePreTrans(x, eig_min);
500 }
501
502 void
solvePre(Vector & x,double & eig_min)503 QuasiTriangular::solvePre(Vector &x, double &eig_min)
504 {
505 addUnit();
506 for (diag_iter di = diag_begin(); di != diag_end(); ++di)
507 {
508 double eig_size;
509 if (!di->isReal())
510 {
511 eig_size = di->getDeterminant();
512 eliminateLeft(di->getIndex()+1, di->getIndex(), x);
513 }
514 else
515 eig_size = *di->getAlpha()*(*di->getAlpha());
516 eig_min = std::min(eig_min, eig_size);
517 }
518
519 blas_int nn = diagonal.getSize();
520 blas_int lda = ld;
521 blas_int incx = x.skip();
522 dtrsv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
523 }
524
525 void
solvePreTrans(Vector & x,double & eig_min)526 QuasiTriangular::solvePreTrans(Vector &x, double &eig_min)
527 {
528 addUnit();
529 for (diag_iter di = diag_begin(); di != diag_end(); ++di)
530 {
531 double eig_size;
532 if (!di->isReal())
533 {
534 eig_size = di->getDeterminant();
535 eliminateRight(di->getIndex()+1, di->getIndex(), x);
536 }
537 else
538 eig_size = *di->getAlpha()*(*di->getAlpha());
539 if (eig_size < eig_min)
540 eig_min = eig_size;
541 }
542
543 blas_int nn = diagonal.getSize();
544 blas_int lda = ld;
545 blas_int incx = x.skip();
546 dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
547 }
548
549 // Calculates x = T·b
550 void
multVec(Vector & x,const ConstVector & b) const551 QuasiTriangular::multVec(Vector &x, const ConstVector &b) const
552 {
553 x = b;
554 blas_int nn = diagonal.getSize();
555 blas_int lda = ld;
556 blas_int incx = x.skip();
557 dtrmv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
558 for (const_diag_iter di = diag_begin(); di != diag_end(); ++di)
559 if (!di->isReal())
560 {
561 int jbar = di->getIndex();
562 x[jbar+1] += di->getBeta2()*(b[jbar]);
563 }
564 }
565
566 void
multVecTrans(Vector & x,const ConstVector & b) const567 QuasiTriangular::multVecTrans(Vector &x, const ConstVector &b) const
568 {
569 x = b;
570 blas_int nn = diagonal.getSize();
571 blas_int lda = ld;
572 blas_int incx = x.skip();
573 dtrmv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
574 for (const_diag_iter di = diag_begin(); di != diag_end(); ++di)
575 if (!di->isReal())
576 {
577 int jbar = di->getIndex();
578 x[jbar] += di->getBeta2()*b[jbar+1];
579 }
580 }
581
582 void
multaVec(Vector & x,const ConstVector & b) const583 QuasiTriangular::multaVec(Vector &x, const ConstVector &b) const
584 {
585 Vector tmp(const_cast<const Vector &>(x)); // new copy
586 multVec(x, b);
587 x.add(1.0, tmp);
588 }
589
590 void
multaVecTrans(Vector & x,const ConstVector & b) const591 QuasiTriangular::multaVecTrans(Vector &x, const ConstVector &b) const
592 {
593 Vector tmp(const_cast<const Vector &>(x)); // new copy
594 multVecTrans(x, b);
595 x.add(1.0, tmp);
596 }
597
598 // Calculates x=x+(this⊗I)·b, where size of I is given by b (KronVector)
599 void
multaKron(KronVector & x,const ConstKronVector & b) const600 QuasiTriangular::multaKron(KronVector &x, const ConstKronVector &b) const
601 {
602 int id = b.getN()*power(b.getM(), b.getDepth()-1);
603 ConstGeneralMatrix b_resh(b, id, b.getM());
604 GeneralMatrix x_resh(x, id, b.getM());
605 x_resh.multAndAdd(b_resh, ConstGeneralMatrix(*this), "trans");
606 }
607
608 // Calculates x=x+(this⊗I)·b, where size of I is given by b (KronVector)
609 void
multaKronTrans(KronVector & x,const ConstKronVector & b) const610 QuasiTriangular::multaKronTrans(KronVector &x, const ConstKronVector &b) const
611 {
612 int id = b.getN()*power(b.getM(), b.getDepth()-1);
613 ConstGeneralMatrix b_resh(b, id, b.getM());
614 GeneralMatrix x_resh(x, id, b.getM());
615 x_resh.multAndAdd(b_resh, ConstGeneralMatrix(*this));
616 }
617
618 void
multKron(KronVector & x) const619 QuasiTriangular::multKron(KronVector &x) const
620 {
621 KronVector b(const_cast<const KronVector &>(x)); // make copy
622 x.zeros();
623 multaKron(x, b);
624 }
625
626 void
multKronTrans(KronVector & x) const627 QuasiTriangular::multKronTrans(KronVector &x) const
628 {
629 KronVector b(const_cast<const KronVector &>(x)); // make copy
630 x.zeros();
631 multaKronTrans(x, b);
632 }
633
634 void
multLeftOther(GeneralMatrix & a) const635 QuasiTriangular::multLeftOther(GeneralMatrix &a) const
636 {
637 a.multLeft(*this);
638 }
639
640 void
multLeftOtherTrans(GeneralMatrix & a) const641 QuasiTriangular::multLeftOtherTrans(GeneralMatrix &a) const
642 {
643 a.multLeftTrans(*this);
644 }
645
646 void
swapDiagLogically(diag_iter it)647 QuasiTriangular::swapDiagLogically(diag_iter it)
648 {
649 diagonal.swapLogically(it);
650 }
651
652 void
checkDiagConsistency(diag_iter it)653 QuasiTriangular::checkDiagConsistency(diag_iter it)
654 {
655 diagonal.checkConsistency(it);
656 }
657
658 double
getAverageDiagSize(diag_iter start,diag_iter end)659 QuasiTriangular::getAverageDiagSize(diag_iter start, diag_iter end)
660 {
661 return diagonal.getAverageSize(start, end);
662 }
663
664 QuasiTriangular::diag_iter
findClosestDiagBlock(diag_iter start,diag_iter end,double a)665 QuasiTriangular::findClosestDiagBlock(diag_iter start, diag_iter end, double a)
666 {
667 return diagonal.findClosestBlock(start, end, a);
668 }
669
670 QuasiTriangular::diag_iter
findNextLargerBlock(diag_iter start,diag_iter end,double a)671 QuasiTriangular::findNextLargerBlock(diag_iter start, diag_iter end, double a)
672 {
673 return diagonal.findNextLargerBlock(start, end, a);
674 }
675
676 int
getNumOffdiagonal() const677 QuasiTriangular::getNumOffdiagonal() const
678 {
679 return diagonal.getSize()*(diagonal.getSize()-1)/2 - diagonal.getNumComplex();
680 }
681