// matrix/sp-matrix.cc
  
  // Copyright 2009-2011  Lukas Burget;  Ondrej Glembek;  Microsoft Corporation
  //                      Saarland University;   Petr Schwarz;   Yanmin Qian;
  //                      Haihua Xu
  
  // See ../../COPYING for clarification regarding multiple authors
  //
  // Licensed under the Apache License, Version 2.0 (the "License");
  // you may not use this file except in compliance with the License.
  // You may obtain a copy of the License at
  
  //  http://www.apache.org/licenses/LICENSE-2.0
  
  // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
  // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
  // MERCHANTABLITY OR NON-INFRINGEMENT.
  // See the Apache 2 License for the specific language governing permissions and
  // limitations under the License.
  
  #include <limits>
  
  #include "matrix/sp-matrix.h"
  #include "matrix/kaldi-vector.h"
  #include "matrix/kaldi-matrix.h"
  #include "matrix/matrix-functions.h"
  #include "matrix/cblas-wrappers.h"
  
  namespace kaldi {
  
  // ****************************************************************************
  // Returns the log-determinant if +ve definite, else KALDI_ERR.
  // ****************************************************************************
  template<typename Real>
  Real SpMatrix<Real>::LogPosDefDet() const {
    TpMatrix<Real> chol(this->NumRows());
    double det = 0.0;
    double diag;
    chol.Cholesky(*this);  // Will throw exception if not +ve definite!
  
    for (MatrixIndexT i = 0; i < this->NumRows(); i++) {
      diag = static_cast<double>(chol(i, i));
      det += kaldi::Log(diag);
    }
    return static_cast<Real>(2*det);
  }
  
  
  template<typename Real>
  void SpMatrix<Real>::Swap(SpMatrix<Real> *other) {
    std::swap(this->data_, other->data_);
    std::swap(this->num_rows_, other->num_rows_);
  }
  
  template<typename Real>
  void SpMatrix<Real>::SymPosSemiDefEig(VectorBase<Real> *s,
                                        MatrixBase<Real> *P,
                                        Real tolerance) const {
    Eig(s, P);
    Real max = s->Max(), min = s->Min();
    KALDI_ASSERT(-min <= tolerance * max);
    s->ApplyFloor(0.0);
  }
  
  template<typename Real>
  Real SpMatrix<Real>::MaxAbsEig() const {
    Vector<Real> s(this->NumRows());
    this->Eig(&s, static_cast<MatrixBase<Real>*>(NULL));
    return std::max(s.Max(), -s.Min());
  }
  
  // returns true if positive definite--uses cholesky.
  template<typename Real>
  bool SpMatrix<Real>::IsPosDef() const {
    MatrixIndexT D = (*this).NumRows();
    KALDI_ASSERT(D > 0);
    try {
      TpMatrix<Real> C(D);
      C.Cholesky(*this);
      for (MatrixIndexT r = 0; r < D; r++)
        if (C(r, r) == 0.0) return false;
      return true;
    }
    catch(...) {  // not positive semidefinite.
      return false;
    }
  }
  
  template<typename Real>
  void SpMatrix<Real>::ApplyPow(Real power) {
    if (power == 1) return;  // can do nothing.
    MatrixIndexT D = this->NumRows();
    KALDI_ASSERT(D > 0);
    Matrix<Real> U(D, D);
    Vector<Real> l(D);
    (*this).SymPosSemiDefEig(&l, &U);
  
    Vector<Real> l_copy(l);
    try {
      l.ApplyPow(power * 0.5);
    }
    catch(...) {
      KALDI_ERR << "Error taking power " << (power * 0.5) << " of vector "
                << l_copy;
    }
    U.MulColsVec(l);
    (*this).AddMat2(1.0, U, kNoTrans, 0.0);
  }
  
  template<typename Real>
  void SpMatrix<Real>::CopyFromMat(const MatrixBase<Real> &M,
                                   SpCopyType copy_type) {
    KALDI_ASSERT(this->NumRows() == M.NumRows() && M.NumRows() == M.NumCols());
    MatrixIndexT D = this->NumRows();
  
    switch (copy_type) {
      case kTakeMeanAndCheck:
        {
          Real good_sum = 0.0, bad_sum = 0.0;
          for (MatrixIndexT i = 0; i < D; i++) {
            for (MatrixIndexT j = 0; j < i; j++) {
              Real a = M(i, j), b = M(j, i), avg = 0.5*(a+b), diff = 0.5*(a-b);
              (*this)(i, j) = avg;
              good_sum += std::abs(avg);
              bad_sum += std::abs(diff);
            }
            good_sum += std::abs(M(i, i));
            (*this)(i, i) = M(i, i);
          }
          if (bad_sum > 0.01 * good_sum) {
            KALDI_ERR << "SpMatrix::Copy(), source matrix is not symmetric: "
                      << bad_sum <<  ">" << good_sum;
          }
          break;
        }
      case kTakeMean:
        {
          for (MatrixIndexT i = 0; i < D; i++) {
            for (MatrixIndexT j = 0; j < i; j++) {
              (*this)(i, j) = 0.5*(M(i, j) + M(j, i));
            }
            (*this)(i, i) = M(i, i);
          }
          break;
        }
      case kTakeLower:
        { // making this one a bit more efficient.
          const Real *src = M.Data();
          Real *dest = this->data_;
          MatrixIndexT stride = M.Stride();
          for (MatrixIndexT i = 0; i < D; i++) {
            for (MatrixIndexT j = 0; j <= i; j++)
              dest[j] = src[j];
            dest += i + 1;
            src += stride;
          }
        }
        break;
      case kTakeUpper:
        for (MatrixIndexT i = 0; i < D; i++)
          for (MatrixIndexT j = 0; j <= i; j++)
            (*this)(i, j) = M(j, i);
        break;
      default:
        KALDI_ASSERT("Invalid argument to SpMatrix::CopyFromMat");
    }
  }
  
  template<typename Real>
  Real SpMatrix<Real>::Trace() const {
    const Real *data = this->data_;
    MatrixIndexT num_rows = this->num_rows_;
    Real ans = 0.0;
    for (int32 i = 1; i <= num_rows; i++, data += i)
      ans += *data;
    return ans;
  }
  
  // diagonal update, this <-- this + diag(v)
  template<typename Real>
  template<typename OtherReal>
  void  SpMatrix<Real>::AddDiagVec(const Real alpha, const VectorBase<OtherReal> &v) {
    int32 num_rows = this->num_rows_;
    KALDI_ASSERT(num_rows == v.Dim() && num_rows > 0);
    const OtherReal *src = v.Data();
    Real *dst = this->data_;
    if (alpha == 1.0)
      for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
        *dst += *src;
    else
      for (int32 i = 1; i <= num_rows; i++, src++, dst += i)
        *dst += alpha * *src;
  }
  
  // instantiate the template above.
  template
  void SpMatrix<float>::AddDiagVec(const float alpha,
                                   const VectorBase<double> &v);
  
  template
  void SpMatrix<double>::AddDiagVec(const double alpha,
                                    const VectorBase<float> &v);
  
  template
  void SpMatrix<float>::AddDiagVec(const float alpha,
                                   const VectorBase<float> &v);
  
  template
  void SpMatrix<double>::AddDiagVec(const double alpha,
                                    const VectorBase<double> &v);
  
  template<>
  template<>
  void SpMatrix<double>::AddVec2(const double alpha, const VectorBase<double> &v);
  
  #ifndef HAVE_ATLAS
  template<typename Real>
  void SpMatrix<Real>::Invert(Real *logdet, Real *det_sign, bool need_inverse) {
    // these are CLAPACK types
    KaldiBlasInt   result;
    KaldiBlasInt   rows = static_cast<int>(this->num_rows_);
    KaldiBlasInt*  p_ipiv = new KaldiBlasInt[rows];
    Real *p_work;  // workspace for the lapack function
    void *temp;
    if ((p_work = static_cast<Real*>(
            KALDI_MEMALIGN(16, sizeof(Real) * rows, &temp))) == NULL) {
      delete[] p_ipiv;
      throw std::bad_alloc();
    }
  #ifdef HAVE_OPENBLAS
    memset(p_work, 0, sizeof(Real) * rows); // gets rid of a probably
    // spurious Valgrind warning about jumps depending upon uninitialized values.
  #endif
  
  
    // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
    // of the data. We have a row-wise storage, therefore, we need to "invert"
    clapack_Xsptrf(&rows, this->data_, p_ipiv, &result);
  
  
    KALDI_ASSERT(result >= 0 && "Call to CLAPACK ssptrf_ called with wrong arguments");
  
    if (result > 0) {  // Singular...
      if (det_sign) *det_sign = 0;
      if (logdet) *logdet = -std::numeric_limits<Real>::infinity();
      if (need_inverse) KALDI_ERR << "CLAPACK stptrf_ : factorization failed";
    } else {  // Not singular.. compute log-determinant if needed.
      if (logdet != NULL || det_sign != NULL) {
        Real prod = 1.0, log_prod = 0.0;
        int sign = 1;
        for (int i = 0; i < (int)this->num_rows_; i++) {
          if (p_ipiv[i] > 0) {  // not a 2x2 block...
            // if (p_ipiv[i] != i+1) sign *= -1;  // row swap.
            Real diag = (*this)(i, i);
            prod *= diag;
          } else {  // negative: 2x2 block. [we are in first of the two].
            i++;  // skip over the first of the pair.
            // each 2x2 block...
            Real diag1 = (*this)(i, i), diag2 = (*this)(i-1, i-1),
                offdiag = (*this)(i, i-1);
            Real thisdet = diag1*diag2 - offdiag*offdiag;
            // thisdet == determinant of 2x2 block.
            // The following line is more complex than it looks: there are 2 offsets of
            // 1 that cancel.
            prod *= thisdet;
          }
          if (i == (int)(this->num_rows_-1) || fabs(prod) < 1.0e-10 || fabs(prod) > 1.0e+10) {
            if (prod < 0) { prod = -prod; sign *= -1; }
            log_prod += kaldi::Log(std::abs(prod));
            prod = 1.0;
          }
        }
        if (logdet != NULL) *logdet = log_prod;
        if (det_sign != NULL) *det_sign = sign;
      }
    }
    if (!need_inverse) {
      delete [] p_ipiv;
      KALDI_MEMALIGN_FREE(p_work);
      return;  // Don't need what is computed next.
    }
    // NOTE: Even though "U" is for upper, lapack assumes column-wise storage
    // of the data. We have a row-wise storage, therefore, we need to "invert"
    clapack_Xsptri(&rows, this->data_, p_ipiv, p_work, &result);
  
    KALDI_ASSERT(result >=0 &&
                 "Call to CLAPACK ssptri_ called with wrong arguments");
  
    if (result != 0) {
      KALDI_ERR << "CLAPACK ssptrf_ : Matrix is singular";
    }
  
    delete [] p_ipiv;
    KALDI_MEMALIGN_FREE(p_work);
  }
  #else
  // in the ATLAS case, these are not implemented using a library and we back off to something else.
  template<typename Real>
  void SpMatrix<Real>::Invert(Real *logdet, Real *det_sign, bool need_inverse) {
    Matrix<Real> M(this->NumRows(), this->NumCols());
    M.CopyFromSp(*this);
    M.Invert(logdet, det_sign, need_inverse);
    if (need_inverse)
      for (MatrixIndexT i = 0; i < this->NumRows(); i++)
        for (MatrixIndexT j = 0; j <= i; j++)
          (*this)(i, j) = M(i, j);
  }
  #endif
  
  template<typename Real>
  void SpMatrix<Real>::InvertDouble(Real *logdet, Real *det_sign,
                                    bool inverse_needed) {
    SpMatrix<double> dmat(*this);
    double logdet_tmp, det_sign_tmp;
    dmat.Invert(logdet ? &logdet_tmp : NULL,
                det_sign ? &det_sign_tmp : NULL,
                inverse_needed);
    if (logdet) *logdet = logdet_tmp;
    if (det_sign) *det_sign = det_sign_tmp;
    (*this).CopyFromSp(dmat);
  }
  
  
  
  double TraceSpSp(const SpMatrix<double> &A, const SpMatrix<double> &B) {
    KALDI_ASSERT(A.NumRows() == B.NumRows());
    const double *Aptr = A.Data();
    const double *Bptr = B.Data();
    MatrixIndexT R = A.NumRows();
    MatrixIndexT RR = (R * (R + 1)) / 2;
    double all_twice = 2.0 * cblas_Xdot(RR, Aptr, 1, Bptr, 1);
    // "all_twice" contains twice the vector-wise dot-product... this is
    // what we want except the diagonal elements are represented
    // twice.
    double diag_once = 0.0;
    for (MatrixIndexT row_plus_two = 2; row_plus_two <= R + 1; row_plus_two++) {
      diag_once += *Aptr * *Bptr;
      Aptr += row_plus_two;
      Bptr += row_plus_two;
    }
    return all_twice - diag_once;
  }
  
  
  float TraceSpSp(const SpMatrix<float> &A, const SpMatrix<float> &B) {
    KALDI_ASSERT(A.NumRows() == B.NumRows());
    const float *Aptr = A.Data();
    const float *Bptr = B.Data();
    MatrixIndexT R = A.NumRows();
    MatrixIndexT RR = (R * (R + 1)) / 2;
    float all_twice = 2.0 * cblas_Xdot(RR, Aptr, 1, Bptr, 1);
    // "all_twice" contains twice the vector-wise dot-product... this is
    // what we want except the diagonal elements are represented
    // twice.
    float diag_once = 0.0;
    for (MatrixIndexT row_plus_two = 2; row_plus_two <= R + 1; row_plus_two++) {
      diag_once += *Aptr * *Bptr;
      Aptr += row_plus_two;
      Bptr += row_plus_two;
    }
    return all_twice - diag_once;
  }
  
  
  template<typename Real, typename OtherReal>
  Real TraceSpSp(const SpMatrix<Real> &A, const SpMatrix<OtherReal> &B) {
    KALDI_ASSERT(A.NumRows() == B.NumRows());
    Real ans = 0.0;
    const Real *Aptr = A.Data();
    const OtherReal *Bptr = B.Data();
    MatrixIndexT row, col, R = A.NumRows();
    for (row = 0; row < R; row++) {
      for (col = 0; col < row; col++)
        ans += 2.0 * *(Aptr++) * *(Bptr++);
      ans += *(Aptr++) * *(Bptr++);  // Diagonal.
    }
    return ans;
  }
  
  template
  float TraceSpSp<float, double>(const SpMatrix<float> &A, const SpMatrix<double> &B);
  
  template
  double TraceSpSp<double, float>(const SpMatrix<double> &A, const SpMatrix<float> &B);
  
  
  template<typename Real>
  Real TraceSpMat(const SpMatrix<Real> &A, const MatrixBase<Real> &B) {
    KALDI_ASSERT(A.NumRows() == B.NumRows() && A.NumCols() == B.NumCols() &&
                 "KALDI_ERR: TraceSpMat: arguments have mismatched dimension");
    MatrixIndexT R = A.NumRows();
    Real ans = (Real)0.0;
    const Real *Aptr = A.Data(), *Bptr = B.Data();
    MatrixIndexT bStride = B.Stride();
    for (MatrixIndexT r = 0;r < R;r++) {
      for (MatrixIndexT c = 0;c < r;c++) {
        // ans += A(r, c) * (B(r, c) + B(c, r));
        ans += *(Aptr++) * (Bptr[r*bStride + c] + Bptr[c*bStride + r]);
      }
      // ans += A(r, r) * B(r, r);
      ans += *(Aptr++) * Bptr[r*bStride + r];
    }
    return ans;
  }
  
  template
  float TraceSpMat(const SpMatrix<float> &A, const MatrixBase<float> &B);
  
  template
  double TraceSpMat(const SpMatrix<double> &A, const MatrixBase<double> &B);
  
  
  template<typename Real>
  Real TraceMatSpMat(const MatrixBase<Real> &A, MatrixTransposeType transA,
                     const SpMatrix<Real> &B, const MatrixBase<Real> &C,
                     MatrixTransposeType transC) {
    KALDI_ASSERT((transA == kTrans?A.NumCols():A.NumRows()) ==
                 (transC == kTrans?C.NumRows():C.NumCols()) &&
                 (transA == kTrans?A.NumRows():A.NumCols()) == B.NumRows() &&
                 (transC == kTrans?C.NumCols():C.NumRows()) == B.NumRows() &&
                 "TraceMatSpMat: arguments have wrong dimension.");
    Matrix<Real> tmp(B.NumRows(), B.NumRows());
    tmp.AddMatMat(1.0, C, transC, A, transA, 0.0);  // tmp = C * A.
    return TraceSpMat(B, tmp);
  }
  
  template
  float TraceMatSpMat(const MatrixBase<float> &A, MatrixTransposeType transA,
                      const SpMatrix<float> &B, const MatrixBase<float> &C,
                      MatrixTransposeType transC);
  template
  double TraceMatSpMat(const MatrixBase<double> &A, MatrixTransposeType transA,
                       const SpMatrix<double> &B, const MatrixBase<double> &C,
                       MatrixTransposeType transC);
  
  template<typename Real>
  Real TraceMatSpMatSp(const MatrixBase<Real> &A, MatrixTransposeType transA,
                       const SpMatrix<Real> &B, const MatrixBase<Real> &C,
                       MatrixTransposeType transC, const SpMatrix<Real> &D) {
    KALDI_ASSERT((transA == kTrans ?A.NumCols():A.NumRows() == D.NumCols()) &&
                 (transA == kTrans ? A.NumRows():A.NumCols() == B.NumRows()) &&
                 (transC == kTrans ? A.NumCols():A.NumRows() == B.NumCols()) &&
                 (transC == kTrans ? A.NumRows():A.NumCols() == D.NumRows()) &&
                 "KALDI_ERR: TraceMatSpMatSp: arguments have mismatched dimension.");
    // Could perhaps optimize this more depending on dimensions of quantities.
    Matrix<Real> tmpAB(transA == kTrans ? A.NumCols():A.NumRows(), B.NumCols());
    tmpAB.AddMatSp(1.0, A, transA, B, 0.0);
    Matrix<Real> tmpCD(transC == kTrans ? C.NumCols():C.NumRows(), D.NumCols());
    tmpCD.AddMatSp(1.0, C, transC, D, 0.0);
    return TraceMatMat(tmpAB, tmpCD, kNoTrans);
  }
  
  template
  float TraceMatSpMatSp(const MatrixBase<float> &A, MatrixTransposeType transA,
                        const SpMatrix<float> &B, const MatrixBase<float> &C,
                        MatrixTransposeType transC, const SpMatrix<float> &D);
  template
  double TraceMatSpMatSp(const MatrixBase<double> &A, MatrixTransposeType transA,
                         const SpMatrix<double> &B, const MatrixBase<double> &C,
                         MatrixTransposeType transC, const SpMatrix<double> &D);
  
  
  template<typename Real>
  bool SpMatrix<Real>::IsDiagonal(Real cutoff) const {
    MatrixIndexT R = this->NumRows();
    Real bad_sum = 0.0, good_sum = 0.0;
    for (MatrixIndexT i = 0; i < R; i++) {
      for (MatrixIndexT j = 0; j <= i; j++) {
        if (i == j)
          good_sum += std::abs((*this)(i, j));
        else
          bad_sum += std::abs((*this)(i, j));
      }
    }
    return (!(bad_sum > good_sum * cutoff));
  }
  
  template<typename Real>
  bool SpMatrix<Real>::IsUnit(Real cutoff) const {
    MatrixIndexT R = this->NumRows();
    Real max = 0.0;  // max error
    for (MatrixIndexT i = 0; i < R; i++)
      for (MatrixIndexT j = 0; j <= i; j++)
        max = std::max(max, static_cast<Real>(std::abs((*this)(i, j) -
                                                       (i == j ? 1.0 : 0.0))));
    return (max <= cutoff);
  }
  
  template<typename Real>
  bool SpMatrix<Real>::IsTridiagonal(Real cutoff) const {
    MatrixIndexT R = this->NumRows();
    Real max_abs_2diag = 0.0, max_abs_offdiag = 0.0;
    for (MatrixIndexT i = 0; i < R; i++)
      for (MatrixIndexT j = 0; j <= i; j++) {
        if (j+1 < i)
          max_abs_offdiag = std::max(max_abs_offdiag,
                                     std::abs((*this)(i, j)));
        else
          max_abs_2diag = std::max(max_abs_2diag,
                                   std::abs((*this)(i, j)));
      }
    return (max_abs_offdiag <= cutoff * max_abs_2diag);
  }
  
  template<typename Real>
  bool SpMatrix<Real>::IsZero(Real cutoff) const {
    if (this->num_rows_ == 0) return true;
    return (this->Max() <= cutoff && this->Min() >= -cutoff);
  }
  
  template<typename Real>
  Real SpMatrix<Real>::FrobeniusNorm() const {
    Real sum = 0.0;
    MatrixIndexT R = this->NumRows();
    for (MatrixIndexT i = 0; i < R; i++) {
      for (MatrixIndexT j = 0; j < i; j++)
        sum += (*this)(i, j) * (*this)(i, j) * 2;
      sum += (*this)(i, i) * (*this)(i, i);
    }
    return std::sqrt(sum);
  }
  
  template<typename Real>
  bool SpMatrix<Real>::ApproxEqual(const SpMatrix<Real> &other, float tol) const {
    if (this->NumRows() != other.NumRows())
      KALDI_ERR << "SpMatrix::AproxEqual, size mismatch, "
                << this->NumRows() << " vs. " << other.NumRows();
    SpMatrix<Real> tmp(*this);
    tmp.AddSp(-1.0, other);
    return (tmp.FrobeniusNorm() <= tol * std::max(this->FrobeniusNorm(), other.FrobeniusNorm()));
  }
  
  // function Floor: A = Floor(B, alpha * C) ... see tutorial document.
  template<typename Real>
  int SpMatrix<Real>::ApplyFloor(const SpMatrix<Real> &C, Real alpha,
                                 bool verbose) {
    MatrixIndexT dim = this->NumRows();
    int nfloored = 0;
    KALDI_ASSERT(C.NumRows() == dim);
    KALDI_ASSERT(alpha > 0);
    TpMatrix<Real> L(dim);
    L.Cholesky(C);
    L.Scale(std::sqrt(alpha));  // equivalent to scaling C by alpha.
    TpMatrix<Real> LInv(L);
    LInv.Invert();
  
    SpMatrix<Real> D(dim);
    {  // D = L^{-1} * (*this) * L^{-T}
      Matrix<Real> LInvFull(LInv);
      D.AddMat2Sp(1.0, LInvFull, kNoTrans, (*this), 0.0);
    }
  
    Vector<Real> l(dim);
    Matrix<Real> U(dim, dim);
  
    D.Eig(&l, &U);
  
    if (verbose) {
      KALDI_LOG << "ApplyFloor: flooring following diagonal to 1: " << l;
    }
    for (MatrixIndexT i = 0; i < l.Dim(); i++) {
      if (l(i) < 1.0) {
        nfloored++;
        l(i) = 1.0;
      }
    }
    l.ApplyPow(0.5);
    U.MulColsVec(l);
    D.AddMat2(1.0, U, kNoTrans, 0.0);
    {  // D' := U * diag(l') * U^T ... l'=floor(l, 1)
      Matrix<Real> LFull(L);
      (*this).AddMat2Sp(1.0, LFull, kNoTrans, D, 0.0);  // A := L * D' * L^T
    }
    return nfloored;
  }
  
  template<typename Real>
  Real SpMatrix<Real>::LogDet(Real *det_sign) const {
    Real log_det;
    SpMatrix<Real> tmp(*this);
    // false== output not needed (saves some computation).
    tmp.Invert(&log_det, det_sign, false);
    return log_det;
  }
  
  
  template<typename Real>
  int SpMatrix<Real>::ApplyFloor(Real floor) {
    MatrixIndexT Dim = this->NumRows();
    int nfloored = 0;
    Vector<Real> s(Dim);
    Matrix<Real> P(Dim, Dim);
    (*this).Eig(&s, &P);
    for (MatrixIndexT i = 0; i < Dim; i++) {
      if (s(i) < floor) {
        nfloored++;
        s(i) = floor;
      }
    }
    (*this).AddMat2Vec(1.0, P, kNoTrans, s, 0.0);
    return nfloored;
  }
  
  template<typename Real>
  MatrixIndexT SpMatrix<Real>::LimitCond(Real maxCond, bool invert) {  // e.g. maxCond = 1.0e+05.
    MatrixIndexT Dim = this->NumRows();
    Vector<Real> s(Dim);
    Matrix<Real> P(Dim, Dim);
    (*this).SymPosSemiDefEig(&s, &P);
    KALDI_ASSERT(maxCond > 1);
    Real floor = s.Max() / maxCond;
    if (floor < 0) floor = 0;
    if (floor < 1.0e-40) {
      KALDI_WARN << "LimitCond: limiting " << floor << " to 1.0e-40";
      floor = 1.0e-40;
    }
    MatrixIndexT nfloored = 0;
    for (MatrixIndexT i = 0; i < Dim; i++) {
      if (s(i) <= floor) nfloored++;
      if (invert)
        s(i) = 1.0 / std::sqrt(std::max(s(i), floor));
      else
        s(i) = std::sqrt(std::max(s(i), floor));
    }
    P.MulColsVec(s);
    (*this).AddMat2(1.0, P, kNoTrans, 0.0);  // (*this) = P*P^T.  ... (*this) = P * floor(s) * P^T  ... if P was original P.
    return nfloored;
  }
  
  void SolverOptions::Check() const {
    KALDI_ASSERT(K>10 && eps<1.0e-10);
  }
  
  template<> double SolveQuadraticProblem(const SpMatrix<double> &H,
                                          const VectorBase<double> &g,
                                          const SolverOptions &opts,
                                          VectorBase<double> *x) {
    KALDI_ASSERT(H.NumRows() == g.Dim() && g.Dim() == x->Dim() && x->Dim() != 0);
    opts.Check();
    MatrixIndexT dim = x->Dim();
    if (H.IsZero(0.0)) {
      KALDI_WARN << "Zero quadratic term in quadratic vector problem for "
                 << opts.name << ": leaving it unchanged.";
      return 0.0;
    }
    if (opts.diagonal_precondition) {
      // We can re-cast the problem with a diagonal preconditioner to
      // make H better-conditioned.
      Vector<double> H_diag(dim);
      H_diag.CopyDiagFromSp(H);
      H_diag.ApplyFloor(std::numeric_limits<double>::min() * 1.0E+3);
      Vector<double> H_diag_sqrt(H_diag);
      H_diag_sqrt.ApplyPow(0.5);
      Vector<double> H_diag_inv_sqrt(H_diag_sqrt);
      H_diag_inv_sqrt.InvertElements();
      Vector<double> x_scaled(*x);
      x_scaled.MulElements(H_diag_sqrt);
      Vector<double> g_scaled(g);
      g_scaled.MulElements(H_diag_inv_sqrt);
      SpMatrix<double> H_scaled(dim);
      H_scaled.AddVec2Sp(1.0, H_diag_inv_sqrt, H, 0.0);
      double ans;
      SolverOptions new_opts(opts);
      new_opts.diagonal_precondition = false;
      ans = SolveQuadraticProblem(H_scaled, g_scaled, new_opts, &x_scaled);
      x->CopyFromVec(x_scaled);
      x->MulElements(H_diag_inv_sqrt);
      return ans;
    }
    Vector<double> gbar(g);
    if (opts.optimize_delta) gbar.AddSpVec(-1.0, H, *x, 1.0);  // gbar = g - H x
    Matrix<double> U(dim, dim);
    Vector<double> l(dim);
    H.SymPosSemiDefEig(&l, &U);  // does svd H = U L V^T and checks that H == U L U^T to within a tolerance.
    // floor l.
    double f = std::max(static_cast<double>(opts.eps), l.Max() / opts.K);
    MatrixIndexT nfloored = 0;
    for (MatrixIndexT i = 0; i < dim; i++) {  // floor l.
      if (l(i) < f) {
        nfloored++;
        l(i) = f;
      }
    }
    if (nfloored != 0 && opts.print_debug_output) {
      KALDI_LOG << "Solving quadratic problem for " << opts.name
                << ": floored " << nfloored<< " eigenvalues. ";
    }
    Vector<double> tmp(dim);
    tmp.AddMatVec(1.0, U, kTrans, gbar, 0.0);  // tmp = U^T \bar{g}
    tmp.DivElements(l);  // divide each element of tmp by l: tmp = \tilde{L}^{-1} U^T \bar{g}
    Vector<double> delta(dim);
    delta.AddMatVec(1.0, U, kNoTrans, tmp, 0.0);  // delta = U tmp = U \tilde{L}^{-1} U^T \bar{g}
    Vector<double> &xhat(tmp);
    xhat.CopyFromVec(delta);
    if (opts.optimize_delta) xhat.AddVec(1.0, *x);  // xhat = x + delta.
    double auxf_before = VecVec(g, *x) - 0.5 * VecSpVec(*x, H, *x),
           auxf_after = VecVec(g, xhat) - 0.5 * VecSpVec(xhat, H, xhat);
    if (auxf_after < auxf_before) {  // Reject change.
      if (auxf_after < auxf_before - 1.0e-10 && opts.print_debug_output)
        KALDI_WARN << "Optimizing vector auxiliary function for "
                   << opts.name<< ": auxf decreased " << auxf_before
                   << " to "  << auxf_after <<  ", change is "
                   << (auxf_after-auxf_before);
      return 0.0;
    } else {
      x->CopyFromVec(xhat);
      return auxf_after - auxf_before;
    }
  }
  
  template<> float SolveQuadraticProblem(const SpMatrix<float> &H,
                                         const VectorBase<float> &g,
                                         const SolverOptions &opts,
                                         VectorBase<float> *x) {
    KALDI_ASSERT(H.NumRows() == g.Dim() && g.Dim() == x->Dim() && x->Dim() != 0);
    SpMatrix<double> Hd(H);
    Vector<double> gd(g);
    Vector<double> xd(*x);
    float ans = static_cast<float>(SolveQuadraticProblem(Hd, gd, opts, &xd));
    x->CopyFromVec(xd);
    return ans;
  }
  
  // Maximizes the auxiliary function   Q(x) = tr(M^T SigmaInv Y) - 0.5 tr(SigmaInv M Q M^T).
  // Like a numerically stable version of   M := Y Q^{-1}.
  template<typename Real>
  Real
  SolveQuadraticMatrixProblem(const SpMatrix<Real> &Q,
                              const MatrixBase<Real> &Y,
                              const SpMatrix<Real> &SigmaInv,
                              const SolverOptions &opts,
                              MatrixBase<Real> *M) {
    KALDI_ASSERT(Q.NumRows() == M->NumCols() &&
                 SigmaInv.NumRows() == M->NumRows() && Y.NumRows() == M->NumRows()
                 && Y.NumCols() == M->NumCols() && M->NumCols() != 0);
    opts.Check();
    MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
    if (Q.IsZero(0.0)) {
      KALDI_WARN << "Zero quadratic term in quadratic matrix problem for "
                 << opts.name << ": leaving it unchanged.";
      return 0.0;
    }
  
    if (opts.diagonal_precondition) {
      // We can re-cast the problem with a diagonal preconditioner in the space
      // of Q (columns of M).  Helps to improve the condition of Q.
      Vector<Real> Q_diag(cols);
      Q_diag.CopyDiagFromSp(Q);
      Q_diag.ApplyFloor(std::numeric_limits<Real>::min() * 1.0E+3);
      Vector<Real> Q_diag_sqrt(Q_diag);
      Q_diag_sqrt.ApplyPow(0.5);
      Vector<Real> Q_diag_inv_sqrt(Q_diag_sqrt);
      Q_diag_inv_sqrt.InvertElements();
      Matrix<Real> M_scaled(*M);
      M_scaled.MulColsVec(Q_diag_sqrt);
      Matrix<Real> Y_scaled(Y);
      Y_scaled.MulColsVec(Q_diag_inv_sqrt);
      SpMatrix<Real> Q_scaled(cols);
      Q_scaled.AddVec2Sp(1.0, Q_diag_inv_sqrt, Q, 0.0);
      Real ans;
      SolverOptions new_opts(opts);
      new_opts.diagonal_precondition = false;
      ans = SolveQuadraticMatrixProblem(Q_scaled, Y_scaled, SigmaInv,
                                        new_opts, &M_scaled);
      M->CopyFromMat(M_scaled);
      M->MulColsVec(Q_diag_inv_sqrt);
      return ans;
    }
  
    Matrix<Real> Ybar(Y);
    if (opts.optimize_delta) {
      Matrix<Real> Qfull(Q);
      Ybar.AddMatMat(-1.0, *M, kNoTrans, Qfull, kNoTrans, 1.0);
    } // Ybar = Y - M Q.
    Matrix<Real> U(cols, cols);
    Vector<Real> l(cols);
    Q.SymPosSemiDefEig(&l, &U);  // does svd Q = U L V^T and checks that Q == U L U^T to within a tolerance.
    // floor l.
    Real f = std::max<Real>(static_cast<Real>(opts.eps), l.Max() / opts.K);
    MatrixIndexT nfloored = 0;
    for (MatrixIndexT i = 0; i < cols; i++) {  // floor l.
      if (l(i) < f) { nfloored++; l(i) = f; }
    }
    if (nfloored != 0 && opts.print_debug_output)
      KALDI_LOG << "Solving matrix problem for " << opts.name
                << ": floored " << nfloored << " eigenvalues. ";
    Matrix<Real> tmpDelta(rows, cols);
    tmpDelta.AddMatMat(1.0, Ybar, kNoTrans, U, kNoTrans, 0.0);  // tmpDelta = Ybar * U.
    l.InvertElements(); KALDI_ASSERT(1.0/l.Max() != 0);  // check not infinite.  eps should take care of this.
    tmpDelta.MulColsVec(l);  // tmpDelta = Ybar * U * \tilde{L}^{-1}
  
    Matrix<Real> Delta(rows, cols);
    Delta.AddMatMat(1.0, tmpDelta, kNoTrans, U, kTrans, 0.0);  // Delta = Ybar * U * \tilde{L}^{-1} * U^T
  
    Real auxf_before, auxf_after;
    SpMatrix<Real> MQM(rows);
    Matrix<Real> &SigmaInvY(tmpDelta);
    { Matrix<Real> SigmaInvFull(SigmaInv);  SigmaInvY.AddMatMat(1.0, SigmaInvFull, kNoTrans, Y, kNoTrans, 0.0); }
    {  // get auxf_before.      Q(x) = tr(M^T SigmaInv Y) - 0.5 tr(SigmaInv M Q M^T).
      MQM.AddMat2Sp(1.0, *M, kNoTrans, Q, 0.0);
      auxf_before = TraceMatMat(*M, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
    }
  
    Matrix<Real> Mhat(Delta);
    if (opts.optimize_delta) Mhat.AddMat(1.0, *M);  // Mhat = Delta + M.
  
    {  // get auxf_after.
      MQM.AddMat2Sp(1.0, Mhat, kNoTrans, Q, 0.0);
      auxf_after = TraceMatMat(Mhat, SigmaInvY, kaldi::kTrans) - 0.5*TraceSpSp(SigmaInv, MQM);
    }
  
    if (auxf_after < auxf_before) {
      if (auxf_after < auxf_before - 1.0e-10)
        KALDI_WARN << "Optimizing matrix auxiliary function for "
                   << opts.name << ", auxf decreased "
                   << auxf_before << " to " << auxf_after << ", change is "
                   << (auxf_after-auxf_before);
      return 0.0;
    } else {
      M->CopyFromMat(Mhat);
      return auxf_after - auxf_before;
    }
  }
  
  template<typename Real>
  Real SolveDoubleQuadraticMatrixProblem(const MatrixBase<Real> &G,
                                         const SpMatrix<Real> &P1,
                                         const SpMatrix<Real> &P2,
                                         const SpMatrix<Real> &Q1,
                                         const SpMatrix<Real> &Q2,
                                         const SolverOptions &opts,
                                         MatrixBase<Real> *M) {
    KALDI_ASSERT(Q1.NumRows() == M->NumCols() && P1.NumRows() == M->NumRows() &&
                 G.NumRows() == M->NumRows() && G.NumCols() == M->NumCols() &&
                 M->NumCols() != 0 && Q2.NumRows() == M->NumCols() &&
                 P2.NumRows() == M->NumRows());
    MatrixIndexT rows = M->NumRows(), cols = M->NumCols();
    // The following check should not fail as we stipulate P1, P2 and one of Q1
    // or Q2 must be +ve def and other Q1 or Q2 must be +ve semidef.
    TpMatrix<Real> LInv(rows);
    LInv.Cholesky(P1);
    LInv.Invert();  // Will throw exception if fails.
    SpMatrix<Real> S(rows);
    Matrix<Real> LInvFull(LInv);
    S.AddMat2Sp(1.0, LInvFull, kNoTrans, P2, 0.0);  // S := L^{-1} P_2 L^{-T}
    Matrix<Real> U(rows, rows);
    Vector<Real> d(rows);
    S.SymPosSemiDefEig(&d, &U);
    Matrix<Real> T(rows, rows);
    T.AddMatMat(1.0, U, kTrans, LInvFull, kNoTrans, 0.0);  // T := U^T * L^{-1}
  
  #ifdef KALDI_PARANOID  // checking mainly for errors in the code or math.
    {
      SpMatrix<Real> P1Trans(rows);
      P1Trans.AddMat2Sp(1.0, T, kNoTrans, P1, 0.0);
      KALDI_ASSERT(P1Trans.IsUnit(0.01));
    }
    {
      SpMatrix<Real> P2Trans(rows);
      P2Trans.AddMat2Sp(1.0, T, kNoTrans, P2, 0.0);
      KALDI_ASSERT(P2Trans.IsDiagonal(0.01));
    }
  #endif
  
    Matrix<Real> TInv(T);
    TInv.Invert();
    Matrix<Real> Gdash(rows, cols);
    Gdash.AddMatMat(1.0, T, kNoTrans, G, kNoTrans, 0.0);  // G' = T G
    Matrix<Real> MdashOld(rows, cols);
    MdashOld.AddMatMat(1.0, TInv, kTrans, *M, kNoTrans, 0.0);  // M' = T^{-T} M
    Matrix<Real> MdashNew(MdashOld);
    Real objf_impr = 0.0;
    for (MatrixIndexT n = 0; n < rows; n++) {
      SpMatrix<Real> Qsum(Q1);
      Qsum.AddSp(d(n), Q2);
      SubVector<Real> mdash_n = MdashNew.Row(n);
      SubVector<Real> gdash_n = Gdash.Row(n);
  
      Matrix<Real> QsumInv(Qsum);
      try {
        QsumInv.Invert();
        Real old_objf = VecVec(mdash_n, gdash_n)
            - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
        mdash_n.AddMatVec(1.0, QsumInv, kNoTrans, gdash_n, 0.0); // m'_n := g'_n * (Q_1 + d_n Q_2)^{-1}
        Real new_objf = VecVec(mdash_n, gdash_n)
            - 0.5 * VecSpVec(mdash_n, Qsum, mdash_n);
        if (new_objf < old_objf) {
          if (new_objf < old_objf - 1.0e-05) {
            KALDI_WARN << "In double quadratic matrix problem: objective "
                "function decreasing during optimization of " << opts.name
                << ", " << old_objf << "->" << new_objf << ", change is "
                << (new_objf - old_objf);
            KALDI_ERR << "Auxiliary function decreasing."; // Will be caught.
          } else {  // Reset to old value, didn't improve (very close to optimum).
            MdashNew.Row(n).CopyFromVec(MdashOld.Row(n));
          }
        }
        objf_impr += new_objf - old_objf;
      }
      catch (...) {
        KALDI_WARN << "Matrix inversion or optimization failed during double "
            "quadratic problem, solving for" << opts.name
            << ": trying more stable approach.";
        objf_impr += SolveQuadraticProblem(Qsum, gdash_n, opts, &mdash_n);
      }
    }
    M->AddMatMat(1.0, T, kTrans, MdashNew, kNoTrans, 0.0); // M := T^T M'.
    return objf_impr;
  }
  
  // rank-one update, this <-- this + alpha V V'
  template<>
  template<>
  void SpMatrix<float>::AddVec2(const float alpha, const VectorBase<float> &v) {
    KALDI_ASSERT(v.Dim() == this->NumRows());
    cblas_Xspr(v.Dim(), alpha, v.Data(), 1,
               this->data_);
  }
  
  template<class Real>
  void SpMatrix<Real>::AddVec2Sp(const Real alpha, const VectorBase<Real> &v,
                                 const SpMatrix<Real> &S, const Real beta) {
    KALDI_ASSERT(v.Dim() == this->NumRows() && S.NumRows() == this->NumRows());
    const Real *Sdata = S.Data();
    const Real *vdata = v.Data();
    Real *data = this->data_;
    MatrixIndexT dim = this->num_rows_;
    for (MatrixIndexT r = 0; r < dim; r++)
      for (MatrixIndexT c = 0; c <= r; c++, Sdata++, data++)
        *data = beta * *data + alpha * vdata[r] * vdata[c] * *Sdata;
  }
  
  
  // rank-one update, this <-- this + alpha V V'
  template<>
  template<>
  void SpMatrix<double>::AddVec2(const double alpha, const VectorBase<double> &v) {
    KALDI_ASSERT(v.Dim() == num_rows_);
    cblas_Xspr(v.Dim(), alpha, v.Data(), 1, data_);
  }
  
  
  template<typename Real>
  template<typename OtherReal>
  void SpMatrix<Real>::AddVec2(const Real alpha, const VectorBase<OtherReal> &v) {
    KALDI_ASSERT(v.Dim() == this->NumRows());
    Real *data = this->data_;
    const OtherReal *v_data = v.Data();
    MatrixIndexT nr = this->num_rows_;
    for (MatrixIndexT i = 0; i < nr; i++)
      for (MatrixIndexT j = 0; j <= i; j++, data++)
        *data += alpha * v_data[i] * v_data[j];
  }
  
  // instantiate the template above.
  template
  void SpMatrix<float>::AddVec2(const float alpha, const VectorBase<double> &v);
  template
  void SpMatrix<double>::AddVec2(const double alpha, const VectorBase<float> &v);
  
  
  template<typename Real>
  Real VecSpVec(const VectorBase<Real> &v1, const SpMatrix<Real> &M,
                const VectorBase<Real> &v2) {
    MatrixIndexT D = M.NumRows();
    KALDI_ASSERT(v1.Dim() == D && v1.Dim() == v2.Dim());
    Vector<Real> tmp_vec(D);
    cblas_Xspmv(D, 1.0, M.Data(), v1.Data(), 1, 0.0, tmp_vec.Data(), 1);
    return VecVec(tmp_vec, v2);
  }
  
  template
  float VecSpVec(const VectorBase<float> &v1, const SpMatrix<float> &M,
                 const VectorBase<float> &v2);
  template
  double VecSpVec(const VectorBase<double> &v1, const SpMatrix<double> &M,
                  const VectorBase<double> &v2);
  
  
  template<typename Real>
  void SpMatrix<Real>::AddMat2Sp(
      const Real alpha, const MatrixBase<Real> &M,
      MatrixTransposeType transM, const SpMatrix<Real> &A, const Real beta) {
    if (transM == kNoTrans) {
      KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
    } else {
      KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
    }
    Vector<Real> tmp_vec(A.NumRows());
    Real *tmp_vec_data = tmp_vec.Data();
    SpMatrix<Real> tmp_A;
    const Real *p_A_data = A.Data();
    Real *p_row_data = this->Data();
    MatrixIndexT M_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows()),
        M_same_dim = (transM == kNoTrans ? M.NumRows() : M.NumCols()),
        M_stride = M.Stride(), dim = this->NumRows();
    KALDI_ASSERT(M_same_dim == dim);
  
    const Real *M_data = M.Data();
  
    if (this->Data() <= A.Data() + A.SizeInBytes() &&
        this->Data() + this->SizeInBytes() >= A.Data()) {
      // Matrices A and *this overlap. Make copy of A
      tmp_A.Resize(A.NumRows());
      tmp_A.CopyFromSp(A);
      p_A_data = tmp_A.Data();
    }
  
    if (transM == kNoTrans) {
      for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
        cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.RowData(r), 1, 0.0, tmp_vec_data, 1);
        cblas_Xgemv(transM, r+1, M_other_dim, alpha, M_data, M_stride,
                    tmp_vec_data, 1, beta, p_row_data, 1);
      }
    } else {
      for (MatrixIndexT r = 0; r < dim; r++, p_row_data += r) {
        cblas_Xspmv(A.NumRows(), 1.0, p_A_data, M.Data() + r, M.Stride(), 0.0, tmp_vec_data, 1);
        cblas_Xgemv(transM, M_other_dim, r+1, alpha, M_data, M_stride,
                    tmp_vec_data, 1, beta, p_row_data, 1);
      }
    }
  }
  
  template<typename Real>
  void SpMatrix<Real>::AddSmat2Sp(
      const Real alpha, const MatrixBase<Real> &M,
      MatrixTransposeType transM, const SpMatrix<Real> &A,
      const Real beta) {
    KALDI_ASSERT((transM == kNoTrans && M.NumCols() == A.NumRows()) ||
                 (transM == kTrans && M.NumRows() == A.NumRows()));
    if (transM == kNoTrans) {
      KALDI_ASSERT(M.NumCols() == A.NumRows() && M.NumRows() == this->num_rows_);
    } else {
      KALDI_ASSERT(M.NumRows() == A.NumRows() && M.NumCols() == this->num_rows_);
    }
    MatrixIndexT Adim = A.NumRows(), dim = this->num_rows_;
  
    Matrix<Real> temp_A(A); // represent A as full matrix.
    Matrix<Real> temp_MA(dim, Adim);
    temp_MA.AddSmatMat(1.0, M, transM, temp_A, kNoTrans, 0.0);
  
    // Next-- we want to do *this = alpha * temp_MA * M^T + beta * *this.
    // To make it sparse vector multiplies, since M is sparse, we'd like
    // to do: for each column c, (*this column c) += temp_MA * (M^T's column c.)
    // [ignoring the alpha and beta here.]
    // It's not convenient to process columns in the symmetric
    // packed format because they don't have a constant stride.  However,
    // we can use the fact that temp_MA * M is symmetric, to just assign
    // each row of *this instead of each column.
    // So the final iteration is:
    // for i = 0... dim-1,
    //   [the i'th row of *this] = beta * [the i'th row of *this] + alpha *
    //                               temp_MA * [the i'th column of M].
    // Of course, we only process the first 0 ... i elements of this row,
    // as that's all that are kept in the symmetric packed format.
  
    Matrix<Real> temp_this(*this);
    Real *data = this->data_;
    const Real *Mdata = M.Data(), *MAdata = temp_MA.Data();
    MatrixIndexT temp_MA_stride = temp_MA.Stride(), Mstride = M.Stride();
  
    if (transM == kNoTrans) {
      // The column of M^T corresponds to the rows of the supplied matrix.
      for (MatrixIndexT i = 0; i < dim; i++, data += i) {
        MatrixIndexT num_rows = i + 1, num_cols = Adim;
        Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
                        temp_MA_stride, Mdata + (i * Mstride), 1, beta, data, 1);
      }
    } else {
      // The column of M^T corresponds to the columns of the supplied matrix.
      for (MatrixIndexT i = 0; i < dim; i++, data += i) {
        MatrixIndexT num_rows = i + 1, num_cols = Adim;
        Xgemv_sparsevec(kNoTrans, num_rows, num_cols, alpha, MAdata,
                        temp_MA_stride, Mdata + i, Mstride, beta, data, 1);
      }
    }
  }
  
  template<typename Real>
  void SpMatrix<Real>::AddMat2Vec(const Real alpha,
                                  const MatrixBase<Real> &M,
                                  MatrixTransposeType transM,
                                  const VectorBase<Real> &v,
                                  const Real beta) {
    this->Scale(beta);
    KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows() &&
                  M.NumCols() == v.Dim()) ||
                 (transM == kTrans && this->NumRows() == M.NumCols() &&
                  M.NumRows() == v.Dim()));
  
    if (transM == kNoTrans) {
      const Real *Mdata = M.Data(), *vdata = v.Data();
      Real *data = this->data_;
      MatrixIndexT dim = this->NumRows(), mcols = M.NumCols(),
          mstride = M.Stride();
      for (MatrixIndexT col = 0; col < mcols; col++, vdata++, Mdata += 1)
        cblas_Xspr(dim, *vdata*alpha, Mdata, mstride, data);
    } else {
      const Real *Mdata = M.Data(), *vdata = v.Data();
      Real *data = this->data_;
      MatrixIndexT dim = this->NumRows(), mrows = M.NumRows(),
          mstride = M.Stride();
      for (MatrixIndexT row = 0; row < mrows; row++, vdata++, Mdata += mstride)
        cblas_Xspr(dim, *vdata*alpha, Mdata, 1, data);
    }
  }
  
  template<typename Real>
  void SpMatrix<Real>::AddMat2(const Real alpha, const MatrixBase<Real> &M,
                               MatrixTransposeType transM, const Real beta)  {
    KALDI_ASSERT((transM == kNoTrans && this->NumRows() == M.NumRows())
                 || (transM == kTrans && this->NumRows() == M.NumCols()));
  
    // Cblas has no function *sprk (i.e. symmetric packed rank-k update), so we
    // use as temporary storage a regular matrix of which we only access its lower
    // triangle
  
    MatrixIndexT this_dim = this->NumRows(),
        m_other_dim = (transM == kNoTrans ? M.NumCols() : M.NumRows());
  
    if (this_dim == 0) return;
    if (alpha == 0.0) {
      if (beta != 1.0) this->Scale(beta);
      return;
    }
  
    Matrix<Real> temp_mat(*this); // wastefully copies upper triangle too, but this
    // doesn't dominate O(N) time.
  
    // This function call is hard-coded to update the lower triangle.
    cblas_Xsyrk(transM, this_dim, m_other_dim, alpha, M.Data(),
                M.Stride(), beta, temp_mat.Data(), temp_mat.Stride());
  
    this->CopyFromMat(temp_mat, kTakeLower);
  }
  
  template<typename Real>
  void SpMatrix<Real>::AddTp2Sp(const Real alpha, const TpMatrix<Real> &T,
                                MatrixTransposeType transM, const SpMatrix<Real> &A,
                                const Real beta) {
    Matrix<Real> Tmat(T);
    AddMat2Sp(alpha, Tmat, transM, A, beta);
  }
  
  template<typename Real>
  void SpMatrix<Real>::AddVecVec(const Real alpha, const VectorBase<Real> &v,
                                 const VectorBase<Real> &w) {
    int32 dim = this->NumRows();
    KALDI_ASSERT(dim == v.Dim() && dim == w.Dim() && dim > 0);
    cblas_Xspr2(dim, alpha, v.Data(), 1, w.Data(), 1, this->data_);
  }
  
  
  template<typename Real>
  void SpMatrix<Real>::AddTp2(const Real alpha, const TpMatrix<Real> &T,
                              MatrixTransposeType transM, const Real beta) {
    Matrix<Real> Tmat(T);
    AddMat2(alpha, Tmat, transM, beta);
  }
  
  
  // Explicit instantiation of the class.
  // This needs to be after the definition of all the class member functions.
  
  template class SpMatrix<float>;
  template class SpMatrix<double>;
  
  
  template<typename Real>
  Real TraceSpSpLower(const SpMatrix<Real> &A, const SpMatrix<Real> &B) {
    MatrixIndexT adim = A.NumRows();
    KALDI_ASSERT(adim == B.NumRows());
    MatrixIndexT dim = (adim*(adim+1))/2;
    return cblas_Xdot(dim, A.Data(), 1, B.Data(), 1);
  }
  // Instantiate the template above.
  template
  double TraceSpSpLower(const SpMatrix<double> &A, const SpMatrix<double> &B);
  template
  float TraceSpSpLower(const SpMatrix<float> &A, const SpMatrix<float> &B);
  
  // Instantiate the template above.
  template float SolveQuadraticMatrixProblem(const SpMatrix<float> &Q,
                                             const MatrixBase<float> &Y,
                                             const SpMatrix<float> &SigmaInv,
                                             const SolverOptions &opts,
                                             MatrixBase<float> *M);
  template double SolveQuadraticMatrixProblem(const SpMatrix<double> &Q,
                                              const MatrixBase<double> &Y,
                                              const SpMatrix<double> &SigmaInv,
                                              const SolverOptions &opts,
                                              MatrixBase<double> *M);
  
  // Instantiate the template above.
  template float SolveDoubleQuadraticMatrixProblem(
      const MatrixBase<float> &G,
      const SpMatrix<float> &P1,
      const SpMatrix<float> &P2,
      const SpMatrix<float> &Q1,
      const SpMatrix<float> &Q2,
      const SolverOptions &opts,
      MatrixBase<float> *M);
  
  template double SolveDoubleQuadraticMatrixProblem(
      const MatrixBase<double> &G,
      const SpMatrix<double> &P1,
      const SpMatrix<double> &P2,
      const SpMatrix<double> &Q1,
      const SpMatrix<double> &Q2,
      const SolverOptions &opts,
      MatrixBase<double> *M);
  
  
  
  } // namespace kaldi