ROL_StdLinearOperatorFactory.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef ROL_STDLINEAROPERATORFACTORY_H
#define ROL_STDLINEAROPERATORFACTORY_H

#include "ROL_StdLinearOperator.hpp"

#include <ctime>
#include <cstdlib>
#include <string>

namespace ROL {

template<class Real> 
class StdLinearOperatorFactory {

  template <typename T> using ROL::Ptr = ROL::Ptr<T>;

  typedef LinearOperator<Real>    OP;
  typedef StdLinearOperator<Real> StdOP;

  typedef std::vector<Real>       vector;

private:


  Teuchos::BLAS<int,Real>     blas_;
  ROL::LAPACK<int,Real>   lapack_;

  // Fill x with uniformly-distributed random values from [lower,upper]
  void randomize( vector &x, Real lower=0.0, Real upper=1.0 ) {
    int N = x.size(); 
    for( int i=0; i<N; ++i ) {
      x[i] = lower+(upper-lower)*static_cast<Real>(rand())/static_cast<Real>(RAND_MAX);
    }
  }

  void diagonal( vector &D, const vector& d ) {
    int N = d.size(); 
    int N2 = N*N;
    D.reserve(N2);
    D.assign(N2,0.0);
    for( int i=0; i<N; ++i ) {
      D[(N+1)*i] = d[i];
    }
  }

  // C = A*B with optional transposes
  void multiply( vector &C, const vector &A, const vector &B, bool transA=false, bool transB=false ) {
    int N2 = A.size();
    int N = (std::round(std::sqrt(N2)));
    bool isSquare = N*N == N2;
    ROL_TEST_FOR_EXCEPTION( !isSquare, std::invalid_argument,
      "Error: vector representation A of matrix must have a square "
      "number of elements."); 
    ROL_TEST_FOR_EXCEPTION( B.size() != N2, std::invalid_argument,
      "Error: vectors A and B must have the same length.");
    ROL_TEST_FOR_EXCEPTION( C.size() != N2, std::invalid_argument,
      "Error: vectors A and C must have the same length.");

    char tra = transA : 'T' : 'N';
    char trb = transB : 'T' : 'N';

    blas_.GEMM(tra,trb,N,N,N,1.0,&A[0],&B[0],N,0.0,N);

  }

  // Orthogonalize the columns of A
  void orthogonalize( vector &A ) {
    int N2 = A.size();
    int N = (std::round(std::sqrt(N2)));
    bool isSquare = N*N == N2;
    ROL_TEST_FOR_EXCEPTION( !isSquare, std::invalid_argument,
      "Error: vector representation of matrix must have a square "
      "number of elements."); 

    vector TAU(N,0.0);
    
    int LDA = N;
    Real LWORK1, LWORK2; 
    int INFO = 0;

    // Query workspace 
    lapack_.GEQRF(N,N,&A[0],LDA,&TAU[0],&LWORK1,-1,&INFO);
    ROL_TEST_FOR_EXCEPTION(INFO,std::logic_error,"LAPACK GEQRF LWORK query failed.");
    
    lapack_.ORGQR(N,N,N,&A[0],LDA,&TAU[0],&LWORK2,-1,&INFO);
    ROL_TEST_FOR_EXCEPTION(INFO,std::logic_error,"LAPACK ORGQR LWORK query failed.");

    const int LWORK = std::max(std::abs(LWORK1),std::abs(LWORK2));

    vector WORK(LWORK);
 
    // Factor the input matrix
    lapack_.GEQRF(N,N,&A[0],LDA,&TAU[0],LWORK,&INFO);
    ROL_TEST_FOR_EXCEPTION(INFO,std::logic_error,"LAPACK GEQRF failed with INFO = " << INFO );
    
    // Overwrite the input matrix with the orthogonal matrix Q
    lapack_.ORGQR(N,N,N,&A[0],LDA,&TAU[0],&WORK[0],LWORK,&INFO);    
    ROL_TEST_FOR_EXCEPTION(INFO,std::logic_error,"LAPACK ORGQR failed with INFO = " << INFO );

  }  


public:

  enum class EMatrixType unsigned {
    MATRIX_SPD = 0,    // Symmetric Positive Definite matrix
    MATRIX_SYMMETRIC,  // Symmetric Indefinite matrix
    MATRIX_UNITARY,    // Unitary matrix
    MATRIX_SINGULAR_S, // Singular symmetric matrix
    MATRIX_SINGULAR_N, // Singular nonsymmetric matrix
    MATRIX_DEFAULT,    // Random nonsymmetric matrix
    MATRIX_LAST 
  }; 

  ~StdLinearOperatorFactor(void) {}

  EMatrixType StringToEMatrixType( satd::string s ) {
    s = removeStringFormat(s);
    for ( EMatrixType mt = MATRIX_SPD; mt < MATRIX_LAST; ++ ) {
      if ( !s.compare(removeStringFormat(EStepToString(mt))) ) {
        return mt;
      }
  }

  std::string EMatrixTypeToString( EMatrixType mt ) {
    std::string retString;
    switch( mt ) {
      case MATRIX_SPD:        retString = "Symmetric Positive Definite";  break;
      case MATRIX_SYMMETRIC:  retString = "Symmetric Indefinite";         break;
      case MATRIX_UNITARY:    retString = "Unitary";                      break;
      case MATRIX_SINGULAR_S: retString = "Singular Symmetric";           break;
      case MATRIX_SINGILAR_N: retString = "Singular Nonsymmetric";        break;
      case MATRIX_DEFAULT:    retString = "Default";                      break;
      case MATRIX_LAST:       retString = "Last (dummy)";                 break; 
    }
    return retString;    
  }

  ROL::Ptr<LinearOperator<Real> > getOperator( int size, const std::string &type="" ) const {
    EMatrixType emt = StringToEMatrixType(type);

    

    int n2 = size*size;

    ROL::Ptr<vector> Ap = ROL::makePtr<vector>(n2);

    switch( emt ) {
      case MATRIX_SPD: {
        
        vector d(size);
        randomize(d,1.0,2.0);

        // A = D
        diagonal(*Ap,d);
        
        vector Q(n2);
        randomize(Q,-1.0,1.0);        
        orthogonalize(Q);
 
        // A = D*Q
        multiply(*Ap,*Ap,Q);
        
        // A = Q'*D*Q
        multiply(*Ap,Q,*Ap,true);

      }
      break;

      case MATRIX_SYMMETRIC: {
        vector d(size);
        randomize(d);

        // A = D
        diagonal(*Ap,d);
        
        vector Q(n2);
        randomize(Q,-1.0,1.0);        
        orthogonalize(Q);
 
        // A = D*Q
        multiply(*Ap,*Ap,Q);
        
        // A = Q'*D*Q
        multiply(*Ap,Q,*Ap,true);
      }
      break; 

      case MATRIX_UNITARY: {
        randomize(*Ap);
        orthogonalize(*Ap);  
      } 
  
      case MATRIX_SINGULAR_S: {
        vector d(size);
        randomize(d);
  
        d[0] = 0;

        // A = D
        diagonal(*Ap,d);
        
        vector Q(n2);
        randomize(Q,-1.0,1.0);        
        orthogonalize(Q);
 
        // A = D*Q
        multiply(*Ap,*Ap,Q);
        
        // A = Q'*D*Q
        multiply(*Ap,Q,*Ap,true);
      

      case MATRIX_SINGULAR_N: {

        vector d(size);
        randomize(d,0.0,1.0);
  
        d[0] = 0;

        // A = D
        diagonal(*Ap,d);

        vector V(n2);
        randomize(V,-1.0,1.0);        
        orthogonalize(V);

        // A = D*V'
        multiply(*Ap,*Ap,Q,false,true);
       
        vector U(n2);
        randomize(U,-1.0,1.0);        
        orthogonalize(U);

        // A = U*D*V'
        multiply(*Ap,U,*Ap);
 
      }

      case MATRIX_DEFAULT: 
      default: {
        randomize(*Ap);
      }   

    }
    return ROL::makePtr<StdOP>(Ap);
  }

  

}; // class StdLinearOperatorFactory

} // namespace ROL


#endif // ROL_STDLINEAROPERATORFACTORY_H