la.h

Go to the documentation of this file.

/* 
 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
 * and Daniel Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or
 * modify under the terms of the GNU General Public License as
 * published by Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.  See the text files
 * COPYING and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 *  scythestat/la.h
 *
 */

#ifndef SCYTHE_LA_H
#define SCYTHE_LA_H

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "algorithm.h"
#include "error.h"
#ifdef SCYTHE_LAPACK
#include "lapack.h"
#endif
#else
#include "scythestat/matrix.h"
#include "scythestat/algorithm.h"
#include "scythestat/error.h"
#ifdef SCYTHE_LAPACK
#include "scythestat/lapack.h"
#endif
#endif

#include <numeric>
#include <algorithm>
#include <set>

namespace scythe {

  namespace {
    typedef unsigned int uint;
  }

  /* Matrix transposition */

  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  t (const Matrix<T,PO,PS>& M)
  {
    uint rows = M.rows();
    uint cols = M.cols();
    Matrix<T,RO,Concrete> ret(cols, rows, false);
    if (PO == Col)
      copy<Col,Row>(M, ret);
    else
      copy<Row,Col>(M, ret);

    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  t (const Matrix<T,O,S>& M)
  {
    return t<O,Concrete>(M);
  }

  /* Ones matrix generation */
  
  template <typename T, matrix_order O, matrix_style S>
  Matrix<T,O,S>
  ones (unsigned int rows, unsigned int cols)
  {
    return Matrix<T,O,S> (rows, cols, true, (T) 1);
  }

  template <typename T, matrix_order O>
  Matrix<T, O, Concrete>
  ones (unsigned int rows, unsigned int cols)
  {
    return ones<T, O, Concrete>(rows, cols);
  }

  template <typename T>
  Matrix<T, Col, Concrete>
  ones (unsigned int rows, unsigned int cols)
  {
    return ones<T, Col, Concrete>(rows, cols);
  }

  inline Matrix<double, Col, Concrete>
  ones (unsigned int rows, unsigned int cols)
  {
    return ones<double, Col, Concrete>(rows, cols);
  }

  /* Identity  Matrix generation */

  // This functor contains the working parts of the eye algorithm.
  namespace {
    template <class T> struct eye_alg {
        T operator() (uint i, uint j) {
          if (i == j)
            return (T) 1.0;
          return (T) 0.0;
        }
    };
  }
  
  template <typename T, matrix_order O, matrix_style S>
  Matrix<T,O,S>
  eye (unsigned int k)
  {
    Matrix<T,O,Concrete> ret(k, k, false);
    for_each_ij_set(ret, eye_alg<T>());
    SCYTHE_VIEW_RETURN(T, O, S, ret)
  }

  template <typename T, matrix_order O>
  Matrix<T, O, Concrete>
  eye (uint k)
  {
    return eye<T, O, Concrete>(k);
  }

  template <typename T>
  Matrix<T, Col, Concrete>
  eye (uint k)
  {
    return eye<T, Col, Concrete>(k);
  }

  inline Matrix<double, Col, Concrete>
  eye (uint k)
  {
    return eye<double, Col, Concrete>(k);
  }

  /* Create a k x 1 vector-additive sequence matrix */

  // The seqa algorithm
  namespace {
    template <typename T> struct seqa_alg {
      T cur_; T inc_;
      seqa_alg(T start, T inc) : cur_ (start), inc_ (inc)  {}
      T operator() () { T ret = cur_; cur_ += inc_; return ret; }
    };
  }
  
  template <typename T, matrix_order O, matrix_style S>
  Matrix<T,O,S>
  seqa (T start, T incr, uint rows)
  {
    Matrix<T,O,Concrete> ret(rows, 1, false);
    generate(ret.begin_f(), ret.end_f(), seqa_alg<T>(start, incr));
    SCYTHE_VIEW_RETURN(T, O, S, ret)
  }

  template <typename T, matrix_order O>
  Matrix<T, O, Concrete>
  seqa (T start, T incr, uint rows)
  {
    return seqa<T, O, Concrete>(start, incr, rows);
  }

  template <typename T>
  Matrix<T, Col, Concrete>
  seqa (T start, T incr, uint rows)
  {
    return seqa<T, Col, Concrete>(start, incr, rows);
  }

  inline Matrix<double, Col, Concrete>
  seqa (double start, double incr, uint rows)
  {
    return seqa<double, Col, Concrete>(start, incr, rows);
  }

  /* Uses the STL sort to sort a Matrix in ascending row-major order */
  
  template <matrix_order SORT_ORDER,
            matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  sort (const Matrix<T, PO, PS>& M)
  {
    Matrix<T,RO,Concrete> ret = M;

    std::sort(ret.template begin<SORT_ORDER>(), 
              ret.template end<SORT_ORDER>());

    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

  template <matrix_order SORT_ORDER,
            typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  sort (const Matrix<T,O,S>& M)
  {
    return sort<SORT_ORDER, O, Concrete>(M);
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  sort (const Matrix<T,O,S>& M)
  {
    return sort<O, O, Concrete>(M);
  }

  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  sortc (const Matrix<T, PO, PS>& M)
  {
    Matrix<T,RO,Concrete> ret = M;

    // TODO need to figure out a way to do fully optimized
    // vector iteration
    for (uint col = 0; col < ret.cols(); ++col) {
      Matrix<T,PO,View> column = ret(_, col);
      std::sort(column.begin(), column.end());
    }

    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  sortc(const Matrix<T,O,S>& M)
  {
    return sortc<O,Concrete>(M);
  }

  /* Column bind two matrices */
  
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,RO,RS>
  cbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    SCYTHE_CHECK_10(A.rows() != B.rows(), scythe_conformation_error,
        "Matrices have different numbers of rows");

    Matrix<T,RO,Concrete> ret(A.rows(), A.cols() + B.cols(), false);
    std::copy(B.template begin_f<Col>(), B.template end_f<Col>(),
              std::copy(A.template begin_f<Col>(), 
                        A.template end_f<Col>(), 
                        ret.template begin_f<Col>()));
    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }


  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  cbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    return cbind<PO1,Concrete>(A, B);
  }

  /* Row bind two matrices */
  
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,RO,RS>
  rbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    SCYTHE_CHECK_10(A.cols() != B.cols(), scythe_conformation_error,
        "Matrices have different numbers of columns");

    Matrix<T,RO,Concrete> ret(A.rows() + B.rows(), A.cols(), false);
    std::copy(B.template begin_f<Row>(), B.template end_f<Row>(),
              std::copy(A.template begin_f<Row>(), 
                        A.template end_f<Row>(), 
                        ret.template begin_f<Row>()));
    SCYTHE_VIEW_RETURN(T, RO, RS, ret)
  }

  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  rbind (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B)
  {
    return rbind<PO1,Concrete>(A, B);
  }

  /* Calculates the order of each element in a Matrix */

  // Functor encapsulating the meat of the algorithm
  namespace {
    template <class T,matrix_order O,matrix_style S> struct order_alg {
      Matrix<T,O> M_;
      order_alg (const Matrix<T,O,S>& M) : M_ (M) {}
      uint operator() (T x) {
        Matrix<bool,O> diff = (M_ < x);
        return std::accumulate(diff.begin_f(), diff.end_f(), (uint) 0);
      }
    };
  }
    
  /* NOTE This function used to only work on column vectors.  I see no
   * reason to maintain this restriction.
   */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<unsigned int, RO, RS>
  order (const Matrix<T, PO, PS>& M)
  {
    Matrix<uint, RO, Concrete> ranks(M.rows(), M.cols(), false);
    std::transform(M.begin_f(), M.end_f(), ranks.template begin_f<PO>(),
                   order_alg<T, PO, PS>(M));
    SCYTHE_VIEW_RETURN(uint, RO, RS, ranks)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<unsigned int, O, Concrete>
  order (const Matrix<T,O,S>& M)
  {
    return order<O,Concrete>(M);
  }
  
  /* Selects all the rows of Matrix A for which binary column vector e
   * has an element equal to 1
   */
   
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,RO,RS>
  selif (const Matrix<T,PO1,PS1>& M, const Matrix<bool,PO2,PS2>& e)
  {
    SCYTHE_CHECK_10(M.rows() != e.rows(), scythe_conformation_error,
     "Data matrix and selection vector have different number of rows");
    SCYTHE_CHECK_10(! e.isColVector(), scythe_dimension_error,
        "Selection matrix is not a column vector");

    uint N = std::accumulate(e.begin_f(), e.end_f(), (uint) 0);
    Matrix<T,RO,Concrete> res(N, M.cols(), false);
    int cnt = 0;
    for (uint i = 0; i < e.size(); ++i) {
      if (e[i]) {
        Matrix<T,RO,View> Mvec = M(i, _);
        // TODO again, need optimized vector iteration
        std::copy(Mvec.begin_f(), Mvec.end_f(), 
            res(cnt++, _).begin_f());
      }
    }

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }
 
  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2>
  Matrix<T,PO1,Concrete>
  selif (const Matrix<T,PO1,PS1>& M, const Matrix<bool,PO2,PS2>& e)
  {
    return selif<PO1,Concrete>(M, e);
  }

  /* Find unique elements in a matrix and return a sorted row vector */
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  unique (const Matrix<T, PO, PS>& M)
  {
    std::set<T> u(M.begin_f(), M.end_f());
    Matrix<T,RO,Concrete> res(1, u.size(), false);
    std::copy(u.begin(), u.end(), res.begin_f());

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  unique (const Matrix<T,O,S>& M)
  {
    return unique<O,Concrete>(M);
  }

  /* NOTE I killed reshape.  It seems redundant with resize. DBP */


  /* Make vector out of unique elements of a symmetric Matrix.  
   */
   
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  vech (const Matrix<T, PO, PS>& M)
  {
    SCYTHE_CHECK_30(! M.isSymmetric(), scythe_dimension_error,
        "Matrix not symmetric");

    Matrix<T,RO,Concrete> 
      res((uint) (0.5 * (M.size() - M.rows())) + M.rows(), 1, false);
    typename Matrix<T,RO,Concrete>::forward_iterator it = res.begin_f();

    /* We want to traverse M in storage order if possible so we take
     * the upper triangle of row-order matrices and the lower triangle
     * of column-order matrices.
     */
    if (M.storeorder() == Col) {
      for (uint i = 0; i < M.rows(); ++i) {
        Matrix<T,PO,View> strip = M(i, i, M.rows() - 1, i);
        it = std::copy(strip.begin_f(), strip.end_f(), it);
      }
    } else {
      for (uint j = 0; j < M.cols(); ++j) {
        Matrix<T,PO,View> strip = M(j, j, j, M.cols() - 1);
        it = std::copy(strip.begin_f(), strip.end_f(), it);
      }
    }

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  vech (const Matrix<T,O,S>& M)
  {
    return vech<O,Concrete>(M);
  }

  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  xpnd (const Matrix<T, PO, PS>& v)
  {
    double size_d = -.5 + .5 * std::sqrt(1. + 8 * v.size());
    SCYTHE_CHECK_10(std::fmod(size_d, 1.) != 0., 
        scythe_dimension_error, 
        "Input vector can't generate square matrix");

    uint size = (uint) size_d;
    Matrix<T,RO,Concrete> res(size, size, false);

    /* It doesn't matter if we travel in order here.
     * TODO Might want to use iterators.
     */
    uint cnt = 0;
    for (uint i = 0; i < size; ++i)
      for (uint j = i; j < size; ++j)
        res(i, j) = res(j, i) = v[cnt++];

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  xpnd (const Matrix<T,O,S>& v)
  {
    return xpnd<O,Concrete>(v);
  }

  /* Get the diagonal of a Matrix. */
  
  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  diag (const Matrix<T, PO, PS>& M)
  {
    Matrix<T,RO,Concrete> res(std::min(M.rows(), M.cols()), 1, false);
    
    /* We want to use iterators to maximize speed for both concretes
     * and views, but we always want to tranvers M in order to avoid
     * slowing down concretes.
     */
    uint incr = 1;
    if (PO == Col)
      incr += M.rows();
    else
      incr += M.cols();

    typename Matrix<T,PO,PS>::const_iterator pit;
    typename Matrix<T,RO,Concrete>::forward_iterator rit 
      = res.begin_f();
    for (pit = M.begin(); pit < M.end(); pit += incr)
      *rit++ = *pit;

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  diag (const Matrix<T,O,S>& M)
  {
    return diag<O,Concrete>(M);
  }

  /* Fast calculation of A*B+C. */

  namespace {
    // Algorithm when one matrix is 1x1
    template <matrix_order RO, typename T,
              matrix_order PO1, matrix_style PS1,
              matrix_order PO2, matrix_style PS2>
    void
    gaxpy_alg(Matrix<T,RO,Concrete>& res, const Matrix<T,PO1,PS1>& X,
              const Matrix<T,PO2,PS2>& B, T constant)
    {
      res = Matrix<T,RO,Concrete>(X.rows(), X.cols(), false);
      if (maj_col<RO,PO1,PO2>())
        std::transform(X.template begin_f<Col>(), 
                       X.template end_f<Col>(),
                       B.template begin_f<Col>(),
                       res.template begin_f<Col>(),
                       ax_plus_b<T>(constant));
      else
        std::transform(X.template begin_f<Row>(), 
                       X.template end_f<Row>(),
                       B.template begin_f<Row>(),
                       res.template begin_f<Row>(),
                       ax_plus_b<T>(constant));
    }
  }

  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2,
            matrix_order PO3, matrix_style PS3>
  Matrix<T,RO,RS>
  gaxpy (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B,
         const Matrix<T,PO3,PS3>& C)
  {
    
    Matrix<T, RO, Concrete> res;

    if (A.isScalar() && B.rows() == C.rows() && B.cols() == C.cols()) {
      // Case 1: 1x1 * nXk + nXk
      gaxpy_alg(res, B, C, A[0]);
    } else if (B.isScalar() && A.rows() == C.rows() &&
               A.cols() == C.cols()) {
      // Case 2: m x n  *  1 x 1  +  m x n
      gaxpy_alg(res, A, C, B[0]);
    } else if (A.cols() == B.rows() && A.rows() == C.rows() &&
               B.cols() == C.cols()) {
      // Case 3: m x n  *  n x k  +  m x k

      res = Matrix<T,RO,Concrete> (A.rows(), B.cols(), false);

      /* These are identical to matrix mult, one optimized for
       * row-major and one for col-major.
       */
      
      T tmp;
      if (RO == Col) { // col-major optimized
       for (uint j = 0; j < B.cols(); ++j) {
         for (uint i = 0; i < A.rows(); ++i)
          res(i, j) = C(i, j);
         for (uint l = 0; l < A.cols(); ++l) {
           tmp = B(l, j);
           for (uint i = 0; i < A.rows(); ++i)
             res(i, j) += tmp * A(i, l);
         }
       }
      } else { // row-major optimized
       for (uint i = 0; i < A.rows(); ++i) {
         for (uint j = 0; j < B.cols(); ++j)
           res(i, j) = C(i, j);
         for (uint l = 0; l < B.rows(); ++l) {
           tmp = A(i, l);
           for (uint j = 0; j < B.cols(); ++j)
             res(i, j) += tmp * B(l,j);
         }
       }
      }

    } else {
      SCYTHE_THROW(scythe_conformation_error,
          "Expects (m x n  *  1 x 1  +  m x n)"
          << "or (1 x 1  *  n x k  +  n x k)"
          << "or (m x n  *  n x k  +  m x k)");
    }

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }

  template <typename T, matrix_order PO1, matrix_style PS1,
            matrix_order PO2, matrix_style PS2,
            matrix_order PO3, matrix_style PS3>
  Matrix<T,PO1,Concrete>
  gaxpy (const Matrix<T,PO1,PS1>& A, const Matrix<T,PO2,PS2>& B,
         const Matrix<T,PO3,PS3>& C)
  {
    return gaxpy<PO1,Concrete>(A,B,C);
  }

  template <matrix_order RO, matrix_style RS, typename T,
            matrix_order PO, matrix_style PS>
  Matrix<T, RO, RS>
  crossprod (const Matrix<T, PO, PS>& A)
  {
    /* When rows > 1, we provide differing implementations of the
     * algorithm depending on A's ordering to maximize strided access.
     *
     * The non-vector version of the algorithm fills in a triangle and
     * then copies it over.
     */
    Matrix<T, RO, Concrete> res;
    T tmp;
    if (A.rows() == 1) {
      res = Matrix<T,RO,Concrete>(A.cols(), A.cols(), true);
      for (uint k = 0; k < A.rows(); ++k) {
        for (uint i = 0; i < A.cols(); ++i) {
          tmp = A(k, i);
          for (uint j = i; j < A.cols(); ++j) {
            res(j, i) =
              res(i, j) += tmp * A(k, j);
          }
        }
      }
    } else {
      if (PO == Row) { // row-major optimized
        /* TODO: This is a little slower than the col-major.  Improve.
         */
        res = Matrix<T,RO,Concrete>(A.cols(), A.cols(), true);
        for (uint k = 0; k < A.rows(); ++k) {
          for (uint i = 0; i < A.cols(); ++i) {
            tmp = A(k, i);
            for (uint j = i; j < A.cols(); ++j) {
                res(i, j) += tmp * A(k, j);
            }
          }
        }
        for (uint i = 0; i < A.cols(); ++i)
          for (uint j = i + 1; j < A.cols(); ++j)
            res(j, i) = res(i, j);
      } else { // col-major optimized
        res = Matrix<T,RO,Concrete>(A.cols(), A.cols(), false);
        for (uint j = 0; j < A.cols(); ++j) {
          for (uint i = j; i < A.cols(); ++i) {
            tmp = (T) 0;
            for (uint k = 0; k < A.rows(); ++k)
              tmp += A(k, i) * A(k, j);
            res(i, j) = tmp;
          }
        }
        for (uint i = 0; i < A.cols(); ++i)
          for (uint j = i + 1; j < A.cols(); ++j)
            res(i, j) = res(j, i);
      }
    }

    SCYTHE_VIEW_RETURN(T, RO, RS, res)
  }

  template <typename T, matrix_order O, matrix_style S>
  Matrix<T, O, Concrete>
  crossprod (const Matrix<T,O,S>& M)
  {
    return crossprod<O,Concrete>(M);
  }

#ifdef SCYTHE_LAPACK
  /* Template specializations of for col-major, concrete
   * matrices of doubles that are only available when a lapack library
   * is available.
   */

  template<>
  inline Matrix<>
  gaxpy<Col,Concrete,double,Col,Concrete,Col,Concrete,Col,Concrete>
  (const Matrix<>& A, const Matrix<>& B, const Matrix<>& C)
  {
    SCYTHE_DEBUG_MSG("Using lapack/blas for gaxpy");
    Matrix<> res;

    if (A.isScalar() && B.rows() == C.rows() && B.cols() == C.cols()) {
      // Case 1: 1x1 * nXk + nXk
      gaxpy_alg(res, B, C, A[0]);
    } else if (B.isScalar() && A.rows() == C.rows() &&
               A.cols() == C.cols()) {
      // Case 2: m x n  *  1 x 1  +  m x n
      gaxpy_alg(res, A, C, B[0]);
    } else if (A.cols() == B.rows() && A.rows() == C.rows() &&
               B.cols() == C.cols()) {
      res = C; // NOTE: this copy may eat up speed gains, but can't be
               //       avoided.
      
      // Case 3: m x n  *  n x k  +  m x k
      double* Apnt = A.getArray();
      double* Bpnt = B.getArray();
      double* respnt = res.getArray();
      const double one(1.0);
      int rows = (int) res.rows();
      int cols = (int) res.cols();
      int innerDim = A.cols();

      lapack::dgemm_("N", "N", &rows, &cols, &innerDim, &one, Apnt,
                     &rows, Bpnt, &innerDim, &one, respnt, &rows);

    }
      return res;
  }

  template<>
  inline Matrix<>
  crossprod(const Matrix<>& A)
  {
    SCYTHE_DEBUG_MSG("Using lapack/blas for crossprod");
    // Set up some constants
    const double zero = 0.0;
    const double one = 1.0;

    // Set up return value and arrays
    Matrix<> res(A.cols(), A.cols(), false);
    double* Apnt = A.getArray();
    double* respnt = res.getArray();
    int rows = (int) A.rows();
    int cols = (int) A.cols();

    lapack::dsyrk_("L", "T", &cols, &rows, &one, Apnt, &rows, &zero, respnt,
                   &cols);
    lapack::make_symmetric(respnt, cols);

    return res;
  }

#endif


} // end namespace scythe


#endif /* SCYTHE_LA_H */