Verified Commit f655c170 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Dynare++ / sylvester equation solver: second batch of modernizations

parent 026c0699

dynare++/sylv/cc/BlockDiagonal.cc

+6 −6
Original line number Diff line number Diff line
@@ -4,7 +4,7 @@ 


#include "BlockDiagonal.hh"



#include <cstdio>

#include <iostream>

#include <cstring>



BlockDiagonal::BlockDiagonal(const double *d, int d_size)
@@ -315,7 +315,7 @@ BlockDiagonal::multKronTrans(KronVector &x) const 
void

BlockDiagonal::printInfo() const

{

  printf("Block sizes:");

  std::cout << "Block sizes:";

  int num_blocks = 0;

  const_diag_iter start = diag_begin();

  const_diag_iter end = findBlockStart(start);
@@ -325,14 +325,14 @@ BlockDiagonal::printInfo() const 
      int ei = diagonal.getSize();

      if (end != diag_end())

        ei = (*end).getIndex();

      printf(" %d", ei-si);

      std::cout << ' ' << ei-si;

      num_blocks++;

      start = end;

      end = findBlockStart(start);

    }

  printf("\nNum blocks: %d\n", num_blocks);

  printf("There are %d zeros out of %d\n",

         getNumZeros(), getNumOffdiagonal());

  std::cout << std::endl

            << "Num blocks: " << num_blocks << std::endl

            << "There are " << getNumZeros() << " zeros out of " << getNumOffdiagonal() << std::endl;

}



int

dynare++/sylv/cc/BlockDiagonal.hh

+2 −2
Original line number Diff line number Diff line
@@ -18,12 +18,12 @@ public: 
  BlockDiagonal(int p, const BlockDiagonal &b);

  BlockDiagonal(const BlockDiagonal &b);

  BlockDiagonal(const QuasiTriangular &t);

  const BlockDiagonal &

  BlockDiagonal &

  operator=(const QuasiTriangular &t)

  {

    GeneralMatrix::operator=(t); return *this;

  }

  const BlockDiagonal &operator=(const BlockDiagonal &b);

  BlockDiagonal &operator=(const BlockDiagonal &b);

  ~BlockDiagonal() override

  {

    delete [] row_len; delete [] col_len;

dynare++/sylv/cc/GeneralMatrix.cc

+31 −51
Original line number Diff line number Diff line
@@ -8,11 +8,13 @@ 
#include <dynblas.h>

#include <dynlapack.h>



#include <cstdio>

#include <iostream>

#include <iomanip>

#include <cstring>

#include <cstdlib>

#include <cmath>

#include <limits>

#include <vector>



int GeneralMatrix::md_length = 32;
@@ -80,9 +82,6 @@ GeneralMatrix::GeneralMatrix(const GeneralMatrix &a, const char *dum1, 
  gemm("T", a, "T", b, 1.0, 0.0);

}



GeneralMatrix::~GeneralMatrix()

= default;



void

GeneralMatrix::place(const ConstGeneralMatrix &m, int i, int j)

{
@@ -178,12 +177,10 @@ GeneralMatrix::zeros() 
  if (ld == rows)

    data.zeros();

  else

    {

    for (int i = 0; i < rows; i++)

      for (int j = 0; j < cols; j++)

        get(i, j) = 0;

}

}



void

GeneralMatrix::unit()
@@ -219,12 +216,10 @@ GeneralMatrix::mult(double a) 
  if (ld == rows)

    data.mult(a);

  else

    {

    for (int i = 0; i < rows; i++)

      for (int j = 0; j < cols; j++)

        get(i, j) *= a;

}

}



// here we must be careful for ld

void
@@ -236,12 +231,10 @@ GeneralMatrix::add(double a, const ConstGeneralMatrix &m) 
  if (ld == rows && m.ld == m.rows)

    data.add(a, m.data);

  else

    {

    for (int i = 0; i < rows; i++)

      for (int j = 0; j < cols; j++)

        get(i, j) += a*m.get(i, j);

}

}



void

GeneralMatrix::add(double a, const ConstGeneralMatrix &m, const char *dum)
@@ -396,9 +389,7 @@ void 
ConstGeneralMatrix::multVec(double a, Vector &x, double b, const ConstVector &d) const

{

  if (x.length() != rows || cols != d.length())

    {

    throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");

    }

  if (rows > 0)

    {

      blas_int mm = rows;
@@ -419,9 +410,7 @@ ConstGeneralMatrix::multVecTrans(double a, Vector &x, double b, 
                                 const ConstVector &d) const

{

  if (x.length() != cols || rows != d.length())

    {

    throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");

    }

  if (rows > 0)

    {

      blas_int mm = rows;
@@ -441,24 +430,19 @@ void 
ConstGeneralMatrix::multInvLeft(const char *trans, int mrows, int mcols, int mld, double *d) const

{

  if (rows != cols)

    {

    throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");

    }

  if (cols != mrows)

    {

    throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");

    }



  if (rows > 0)

    {

      GeneralMatrix inv(*this);

      auto *ipiv = new lapack_int[rows];

      std::vector<lapack_int> ipiv(rows);

      lapack_int info;

      lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld;

      dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv, &info);

      dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv, d,

      dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv.data(), &info);

      dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv.data(), d,

             &mld2, &info);

      delete [] ipiv;

    }

}
@@ -481,9 +465,7 @@ void 
ConstGeneralMatrix::multInvLeft(Vector &d) const

{

  if (d.skip() != 1)

    {

    throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");

    }



  multInvLeft("N", d.length(), 1, d.length(), d.base());

}
@@ -493,9 +475,7 @@ void 
ConstGeneralMatrix::multInvLeftTrans(Vector &d) const

{

  if (d.skip() != 1)

    {

    throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");

    }



  multInvLeft("T", d.length(), 1, d.length(), d.base());

}
@@ -523,16 +503,17 @@ ConstGeneralMatrix::isZero() const 
void

ConstGeneralMatrix::print() const

{

  printf("rows=%d, cols=%d\n", rows, cols);

  auto ff = std::cout.flags();

  std::cout << "rows=" << rows << ", cols=" << cols << std::endl;

  for (int i = 0; i < rows; i++)

    {

      printf("row %d:\n", i);

      std::cout << "row " << i << ':' << std::endl

                << std::setprecision(3);

      for (int j = 0; j < cols; j++)

        {

          printf("%6.3g ", get(i, j));

        }

      printf("\n");

        std::cout << std::setw(6) << get(i, j) << ' ';

      std::cout << std::endl;

    }

  std::cout.flags(ff);

}



void
@@ -563,16 +544,15 @@ SVDDecomp::construct(const GeneralMatrix &A) 
  lapack_int lwork = -1;

  lapack_int info;



  auto *iwork = new lapack_int[8*minmn];

  std::vector<lapack_int> iwork(8*minmn);

  // query for optimal lwork

  dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &tmpwork,

         &lwork, iwork, &info);

         &lwork, iwork.data(), &info);

  lwork = (lapack_int) tmpwork;

  Vector work(lwork);

  // do the decomposition

  dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work.base(),

         &lwork, iwork, &info);

  delete [] iwork;

         &lwork, iwork.data(), &info);

  if (info < 0)

    throw SYLV_MES_EXCEPTION("Internal error in SVDDecomp constructor");

  if (info == 0)

dynare++/sylv/cc/GeneralMatrix.hh

+1 −1
Original line number Diff line number Diff line
@@ -145,7 +145,7 @@ public: 
  GeneralMatrix(const GeneralMatrix &a, const char *dum1,

                const GeneralMatrix &b, const char *dum2);



  virtual ~GeneralMatrix();

  virtual ~GeneralMatrix() = default;

  GeneralMatrix &operator=(const GeneralMatrix &m) = default;



  const double &

dynare++/sylv/cc/KronVector.cc

+2 −2
Original line number Diff line number Diff line
@@ -40,7 +40,7 @@ KronVector::KronVector(const ConstKronVector &v) 
  Vector::operator=(v);

}



const KronVector &

KronVector &

KronVector::operator=(const ConstKronVector &v)

{

  Vector::operator=(v);
@@ -50,7 +50,7 @@ KronVector::operator=(const ConstKronVector &v) 
  return *this;

}



const KronVector &

KronVector &

KronVector::operator=(const Vector &v)

{

  if (length() != v.length())