MyLapack.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004 by Claudio Attaccalite                             *
 *   claudio@freescience.info                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
#include "MyLapack.h"

using namespace std;

/********** Double Precision **************/

double DLU_det(DMatrix &A,int *ipiv)
{
        double det=A(0,0);
        int isum=0;
        if(ipiv[0]!=0) isum=1;
        for(int i=1;i<A.dim1();i++)
        {
                det*=A(i,i);
                if(ipiv[i]!=i) isum++;
        }
        if(isum%2!=0) det=-det;
        return det;
}

int DLU_decomp(DMatrix &A,int *ipiv)
{
int info;
int n=A.dim1();
F77_FUNC(dgetrf,DGETRF)(&n,&n,A.array(),&n,ipiv,&info);
return info;
}

int DLU_invert(DMatrix &A,int *ipiv,double *work,int &lwork)
{ 
int info;
int n  =A.dim2();
int lda=A.dim1();
F77_FUNC(dgetri,DGETRI)(&n,A.array(),&lda,ipiv,work,&lwork,&info);
return info;
}


int DSYEV(char jobz,char uplo,DMatrix &M,DVector &V,double *work,int &lwork)
{
int n  =M.dim2();
int lda=M.dim1();
int info;
F77_FUNC(dsyev,DSYEV)(&jobz,&uplo,&n,M.array(),&lda,V.array(),work,&lwork,&info);
 return info;
}

int DLU_solve(DMatrix &M,int *ipiv,DVector &b)
{ 
int info;
int one=1;
int n=M.dim1();
F77_FUNC(dgetrs,DGETRS)("N", &n, &one,M.array(),&n,ipiv,b.array(),&n, &info);
return info;
}

/********** Single Precision **************/

float SLU_det(SMatrix &A,int *ipiv)
{
        float det=A(0,0);
        int isum=0;
        if(ipiv[0]!=0) isum=1;
        for(int i=1;i<A.dim1();i++)
        {
                det*=A(i,i);
                if(ipiv[i]!=i) isum++;
        }
        if(isum%2!=0) det=-det;
        return det;
}

int SLU_decomp(SMatrix &A,int *ipiv)
{
int info;
int n=A.dim1();
F77_FUNC(sgetrf,SGETRF)(&n,&n,A.array(),&n,ipiv,&info);
return info;
}

int SLU_invert(SMatrix &A,int *ipiv,float *work,int &lwork)
{ 
int info;
int n  =A.dim2();
int lda=A.dim1();
F77_FUNC(sgetri,SGETRI)(&n,A.array(),&lda,ipiv,work,&lwork,&info);
return info;
}


int SSYEV(char jobz,char uplo,SMatrix &M,SVector &V,float *work,int &lwork)
{
int n  =M.dim2();
int lda=M.dim1();
int info;
F77_FUNC(ssyev,SSYEV)(&jobz,&uplo,&n,M.array(),&lda,V.array(),work,&lwork,&info);
 return info;
}

int SLU_solve(SMatrix &M,int *ipiv,SVector &b)
{ 
int info;
int one=1;
int n=M.dim1();
F77_FUNC(sgetrs,SGETRS)("N", &n, &one,M.array(),&n,ipiv,b.array(),&n, &info);
return info;
}

---