lapack.h

/*
This file is distributed as part of the BIAS library (Basic ImageAlgorithmS)
but it has not been developed by the authors of BIAS.

For copyright, author and license information see below.
*/


/*
*
* Template Numerical Toolkit (TNT): Linear Algebra Module
*
* Mathematical and Computational Sciences Division
* National Institute of Technology,
* Gaithersburg, MD USA
*
*
* This software was developed at the National Institute of Standards and
* Technology (NIST) by employees of the Federal Government in the course
* of their official duties. Pursuant to title 17 Section 105 of the
* United States Code, this software is not subject to copyright protection
* and is in the public domain.  The Template Numerical Toolkit (TNT) is
* an experimental system.  NIST assumes no responsibility whatsoever for
* its use by other parties, and makes no guarantees, expressed or implied,
* about its quality, reliability, or any other characteristic.
*
* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
* see http://math.nist.gov/tnt for latest updates.
*
*/



// Header file for Fortran Lapack

#ifndef LAPACK_H
#define LAPACK_H

// This file incomplete and included here to only demonstrate the
// basic framework for linking with the Fortran Lapack routines.

#include <Base/Common/W32Compat.hh>
#include "fortran.h"
#include "vec.h"
#include "fmat.h"

using namespace std;


#define F77_DGESV   dgesv_
#define F77_DGELS   dgels_
#define F77_DSYEV   dsyev_

// jw
#define F77_DGEEV   dgeev_
// jw
#define F77_DGESVD  dgesvd_

extern "C"
{

  // linear equations (general) using LU factorizaiton
  //
  void F77_DGESV(cfi_ N, cfi_ nrhs, fda_ A, cfi_ lda,
    fia_ ipiv, fda_ b, cfi_ ldb, fi_ info);

  // solve linear least squares using QR or LU factorization
  //
  void F77_DGELS(cfch_ trans, cfi_ M, 
    cfi_ N, cfi_ nrhs, fda_ A, cfi_ lda, fda_ B, cfi_ ldb, fda_ work, 
    cfi_ lwork, fi_ info);

  // solve symmetric eigenvalues
  //
  void F77_DSYEV( cfch_ jobz, cfch_ uplo, cfi_ N, fda_  A, cfi_ lda, 
    fda_ W, fda_ work, cfi_ lwork, fi_ info);

  // solve unsymmetric eigenvalues
  //
  void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda,
    fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr, 
    cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info);

  /**
  [see fortran.h for typedefs]
  SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO )

  Scalar Arguments:
  CHARACTER          JOBU, JOBVT
  INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N

  Array Arguments:
  DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),  VT( LDVT, * ), WORK( * )

  Purpose
  =======

  DGESVD computes the singular value decomposition (SVD) of a real
  M-by-N matrix A, optionally computing the left and/or right singular
  vectors. The SVD is written

  A = U * SIGMA * transpose(V)

  where SIGMA is an M-by-N matrix which is zero except for its
  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  are the singular values of A; they are real and non-negative, and
  are returned in descending order.  The first min(m,n) columns of
  U and V are the left and right singular vectors of A.

  Note that the routine returns V**T, not V.
  (JW, 03/2002)
  **/
  void F77_DGESVD( cfch_ JOBU,
    cfch_ JOBVT,
    fi_ M, 
    fi_ N, 
    fda_ A,
    fi_ LDA,
    fda_ S,
    fda_ U,
    fi_ LDU,
    fda_ VT,
    fi_ LDVT,
    fda_ WORK,
    fi_ LWORK,
    fi_ INFO  );

}

// solve linear equations using LU factorization

namespace TNT{

  Vector<double> Lapack_LU_linear_solve(const Fortran_Matrix<double> &A,
    const Vector<double> &b)
  {
    const Fortran_integer one=1;
    long int M=A.num_rows();
    long int N=A.num_cols();

    Fortran_Matrix<double> Tmp(A);
    Vector<double> x(b);
    Vector<Fortran_integer> index(M);
    Fortran_integer info = 0;

    F77_DGESV(&N, &one, &Tmp(1,1), &M, &index(1), &x(1), &M, &info);    

    if (info != 0) return Vector<double>(0);
    else
      return x;
  }

  // solve linear least squares problem using QR factorization
  //
  Vector<double> Lapack_LLS_QR_linear_solve(const Fortran_Matrix<double> &A,
    const Vector<double> &b)
  {
    const Fortran_integer one=1;
    Fortran_integer M=A.num_rows();
    Fortran_integer N=A.num_cols();

    Fortran_Matrix<double> Tmp(A);
    Vector<double> x(b);
    Fortran_integer info = 0;

    char transp = 'N';
    Fortran_integer lwork = 5 * (M+N);      // temporary work space
    Vector<double> work(lwork);

    F77_DGELS(&transp, &M, &N, &one, &Tmp(1,1), &M, &x(1), &M,  &work(1),
      &lwork, &info); 

    if (info != 0) return Vector<double>(0);
    else
      return x;
  }

  // *********************** Eigenvalue problems *******************

  // solve symmetric eigenvalue problem (eigenvalues only)
  //
  Vector<double> Upper_symmetric_eigenvalue_solve(const Fortran_Matrix<double> &A)
  {
    char jobz = 'N';
    char uplo = 'U';
    Fortran_integer N = A.num_rows();

    BIASASSERT(N == A.num_cols());

    Vector<double> eigvals(N);
    Fortran_integer worksize = 3*N;
    Fortran_integer info = 0;
    Vector<double> work(worksize);
    Fortran_Matrix<double> Tmp = A;

    F77_DSYEV(&jobz, &uplo, &N, &Tmp(1,1), &N, eigvals.begin(), work.begin(),
      &worksize, &info);

    if (info != 0) return Vector<double>();
    else
      return eigvals;
  }


  // solve unsymmetric eigenvalue problems (for a rectangular matrix)
  //
  int eigenvalue_solve(const Fortran_Matrix<double> &A, 
    Vector<double> &wr, Vector<double> &wi)
  {
    char jobvl = 'N';
    char jobvr = 'N';

    Fortran_integer N = A.num_rows();


    BIASASSERT(N == A.num_cols());

    if (N<1) return 1;

    Fortran_Matrix<double> vl(1,N);  /* should be NxN ? **** */
    Fortran_Matrix<double> vr(1,N);  
    Fortran_integer one = 1;

    Fortran_integer worksize = 5*N;
    Fortran_integer info = 0;
    Vector<double> work(worksize, 0.0);
    Fortran_Matrix<double> Tmp = A;

    wr.newsize(N);
    wi.newsize(N);

    //  void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda,
    //      fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr, 
    //      cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info);

    F77_DGEEV(&jobvl, &jobvr, &N, &Tmp(1,1), &N, &(wr(1)),
      &(wi(1)), &(vl(1,1)), &one, &(vr(1,1)), &one,
      &(work(1)), &worksize, &info);

    return (info==0 ? 0: 1);
  }




  /** solve general eigenproblem for a general quadratix (n x n)  matrix.
  computes eigenvalues and eigenvectors of A.
  (Jan Woetzel), 03/2002
  **/
  int eigenproblem_special_quadratic_matrix_solve(const Fortran_Matrix<double> &A, 
    Vector<double> &ret_EigenValuesRealPart, Vector<double> &ret_EigenValuesImagPart,
    Fortran_Matrix<double> &eigenVecMatrixR ) 
  {        
    char jobvl = 'V'; // calculate the left eigenvectors
    char jobvr = 'V'; // calculate the right eigenvectors

    Fortran_integer N = A.num_rows(); // dim of the Matrix should be NxN
    BIASASSERT(N == A.num_cols() ); // A must be rectangular NxN matrix
    if (N<1) return 1; // 1x1 matrix is a scalar, decomposition is nonsense

    Fortran_integer ldvl = N; // size of the left eigenvectormatrix
    Fortran_integer ldvr = N; // size of the right eigenvectormatrix

    Fortran_Matrix<double> eigenVectsLeft( N,ldvl ); // should be NxN, but left Eigenvectorss are not calculated due to jobvl=N
    Fortran_Matrix<double> eigenVectsRight( N,ldvr ); // will contain a (NxN) Matrix with the right Eigenvectors in columns. 

    Fortran_integer worksize = 5*N;
    Fortran_integer info = 0;
    Vector<double> work(worksize, 0.0);
    Fortran_Matrix<double> Tmp = A;

    ret_EigenValuesRealPart.newsize(N);
    ret_EigenValuesImagPart.newsize(N);

    //  void F77_DGEEV(cfch_ jobvl, cfch_ jobvr, cfi_ N, fda_ A, cfi_ lda,
    //      fda_ wr, fda_ wi, fda_ vl, cfi_ ldvl, fda_ vr, 
    //      cfi_ ldvr, fda_ work, cfi_ lwork, fi_ info);

    // call the lapack lib SVD routine: 
    F77_DGEEV(&jobvl, &jobvr, &N, &Tmp(1,1), &N, &(ret_EigenValuesRealPart(1)),
      &(ret_EigenValuesImagPart(1)), 
      &(eigenVectsLeft(1,1)), &ldvl, 
      &(eigenVectsRight(1,1)), &ldvr,
      &(work(1)), &worksize, &info);


    eigenVecMatrixR = eigenVectsRight; // return a matrix with Eigenvectors in the columns (Spalten) 

    return (info==0 ? 0: 1); // if info !=  0 the the arrays were too small. 
  }




  /** solve general eigenproblem.
  computes singular value decomposition of a rectangular m x n matrix A 
  calls extern liblapack routine. 
  Jan Woetzel 03/2002
  **/
  int singular_value_decomposition_general_rectangular_matrix_solve ( const Fortran_Matrix<double> &A, 
    Vector<double> &ret_S, /* vector with singular values of A with s(i) >= s(i+1) */
    Fortran_Matrix<double> &ret_U, /* m x m orthogonal/unitary matrix U */
    Fortran_Matrix<double> &ret_VT /* n x n orthogonal/unitary matrx transpose(V)*/ ) 
  {   
    char jobu  = 'A'; // calculate all columns of U
    char jobvt = 'A'; // calculate all n rows of V_h

    Fortran_integer M = A.num_rows(); // dim of the Matrix A should be MxN
    Fortran_integer N = A.num_cols(); // dim of the Matrix A should be MxN

    // create a copy od A because jobu=O or jobvt=O may replace the content of A in place !!!
    Fortran_Matrix<double> TmpA = A; // M x N Matrix A to be decomposed 

    Fortran_Matrix<double> U ( M, M );// M x M Orthogonal/Unitary matrix
    Fortran_Matrix<double> V ( N, N );// N x N Orthogonal/Unitary matrix
    Fortran_Matrix<double> VT ( N, N );// N x N Orthogonal/Unitary matrix transpose(V)
    Vector<double> S ( max( (int)M, (int)N) );// the M singular values 

    Fortran_integer lda  = M;
    Fortran_integer ldu  = M; // because jobu = A
    Fortran_integer ldvt = N; // because jobvt = A


    Fortran_integer worksize =  6 * (M > N ? M : N) ;

    Fortran_integer info = 0;
    Vector<double> work(worksize, 0.0);

    // call extern Fortran2C routine from liblapack: 
    F77_DGESVD( &jobu,
      &jobvt,
      &M, 
      &N, 
      &TmpA(1,1),
      &lda,
      &S(1),
      &U(1,1),
      &ldu,
      &VT(1,1),
      &ldvt,
      &work(1),
      &worksize,
      &info);

    //#if 0
    //    // debug output    
    //    cout << " eVals = " << endl;
    //    for (int row=1; row <= M; row++) {
    //  cout <<row<<". : " <<S( row );
    //    };
    //    cout << " U=  " << endl;
    //    for (int row=1; row <= M; row++) { 
    //      for (int col=1; col <= M; col++) { 
    //      cout << U (row, col) << " ";
    //  };
    //  cout<<endl;
    //    };
    //    cout << "  VT=  " << endl;
    //    for (int row=1; row <= N; row++) { 
    //  for (int col=1; col <= N; col++) { 
    //      cout << VT (row, col)<< " " ;
    //  };
    //  cout<<endl;
    //    };
    //#endif

    ret_S  = S; // better compute with the return arguments in place (except A!!) (JW)
    ret_U  = U;
    ret_VT = VT; 

    return (info==0 ? 0: 1); // if info !=  0 the the arrays were too small. 
  }
}


#endif
// LAPACK_H