matrixlapack.c File Reference

wrapper for linear algebra operations using lapack More...

#include "matrixlapack.h"
#include <stdlib.h>
#include "mfmemory.h"

Go to the source code of this file.

Functions
int	dgetrf_ (int *m, int *n, double *a, int *lda, int *ipiv, int *info)
int	dgetri_ (int *n, double *a, int *lda, int *ipiv, double *work, int *lwork, int *info)
int	dsyev_ (char *jobz, char *uplo, int *n, double *a, int *lda, double *w, double *work, int *lwork, int *info)
int	M_DSYEV (char *jobz, char *uplo, int *n, double *a, int *lda, double *w, double *work, int *lwork, int *info)
double	dlamch_ (char *cmach)
int	zheev_ (char *jobz, char *uplo, int *n, double complex *a, int *lda, double *w, double complex *work, int *lwork, double *rwork, int *info)
int	dsyevx_ (char *jobz, char *range, char *uplo, int *n, double *a, int *lda, double *vl, double *vu, int *il, int *iu, double *abstol, int *m, double *w, double *z__, int *ldz, double *work, int *lwork, int *iwork, int *ifail, int *info)
void	to_f (int N, int M, double **A, double *a)
	function for transforming Row-major matrix (C) to Column-major matrix (Fortran) More...
int	DSEVvalue (int xNsize, double **A, double *r)
	obtain eigenvalues of real symmetric A More...
int	DInv (int xNsize, double **xM, double **xIM)
	obtain eigenvalues of inverse of real matrix xM More...
int	DSEVvector (int xNsize, double **A, double *r, double **vec)
	obtain eigenvalues and eigenvectors of real symmetric A More...
int	DSEVXU (int xNsize, double **A, double *r, double **X, int xNev)
	obtain eigenvalues A More...
int	ZHEEVall (int xNsize, double complex **A, double complex *r, double complex **vec)
	obtain eigenvalues and eigenvectors of Hermite matrix A More...

Detailed Description

wrapper for linear algebra operations using lapack

Version

0.1, 0.2

Author

Takahiro Misawa (The University of Tokyo)

Definition in file matrixlapack.c.

Function Documentation

dgetrf_()

int dgetrf_	(	int *	m,
		int *	n,
		double *	a,
		int *	lda,
		int *	ipiv,
		int *	info
	)

Referenced by DInv().

dgetri_()

int dgetri_	(	int *	n,
		double *	a,
		int *	lda,
		int *	ipiv,
		double *	work,
		int *	lwork,
		int *	info
	)

Referenced by DInv().

DInv()

int DInv	(	int	xNsize,
		double **	xM,
		double **	xIM
	)

obtain eigenvalues of inverse of real matrix xM

Parameters

[in]	xNsize	matrix size
[in]	xM	matrix
[out]	xIM	inverse of xM

Returns

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 154 of file matrixlapack.c.

References dgetrf_(), and dgetri_().

                                               {

    int i,j,k;
    int m,n,lda,info,*piv,lwork;
    double *work;
    double *a;

    m=n=lda=lwork=xNsize;

    a = (double*)malloc(xNsize*xNsize*sizeof(double));
    work = (double*)malloc(xNsize*sizeof(double));
    piv = (int*)malloc(xNsize*sizeof(int));

    k=0;
    for(j=0;j<xNsize;j++){
        for(i=0;i<xNsize;i++){
            a[k] = xM[i][j];
            k++;
        }
    }

    dgetrf_(&m, &n, a, &lda, piv, &info);
    dgetri_(&n, a, &lda, piv, work, &lwork, &info);

    if(info != 0){
        free(a);
        free(work);
        free(piv);
        return 0;
    }

    for(k=0;k<xNsize*xNsize;k++){
        xIM[k%xNsize][k/xNsize] = a[k];
    }
    free(a);
    free(work);
    free(piv);

    return 1;
}

dlamch_()

double dlamch_

(

char *

cmach

)

Referenced by DSEVXU().

DSEVvalue()

int DSEVvalue	(	int	xNsize,
		double **	A,
		double *	r
	)

obtain eigenvalues of real symmetric A

Parameters

[in]	xNsize
[in]	A	matrix
[out]	r	eigenvalues

Returns

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 94 of file matrixlapack.c.

References M_DSYEV(), and to_f().

Referenced by Lanczos_EigenValue().

                                                {

  int k;
  char jobz, uplo;
  int n, lda, lwork, info;
  double *a, *w, *work;
#ifdef SR
  int *iwork, liwork;
  liwork = 5 * xNsize + 3;
  iwork = (int*)malloc(liwork * sizeof(double));
#endif

  n = lda = xNsize;
  lwork = 4*xNsize; /* 3*xNsize OK?*/

  a = (double*)malloc(xNsize*xNsize*sizeof(double));
  w = (double*)malloc(xNsize*sizeof(double));
  work = (double*)malloc(lwork*sizeof(double));

  to_f(xNsize, xNsize, A, a);

  jobz = 'N';
  uplo = 'U';

#ifdef SR
  M_DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &iwork, &liwork, &info);
  free(iwork);
#else
  M_DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
#endif

  if(info != 0){
    free(a);
    free(w);
    free(work);
    return 0;
  }

  for(k=0;k<xNsize;k++){
    r[k] = w[k];
  }

  free(a);
  free(w);
  free(work);

  return 1;
}

DSEVvector()

int DSEVvector	(	int	xNsize,
		double **	A,
		double *	r,
		double **	vec
	)

obtain eigenvalues and eigenvectors of real symmetric A

Parameters

xNsize	size of matrix
A	matrix
r	eigenvalues
vec	eignevectos

Returns

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 210 of file matrixlapack.c.

References dsyev_(), M_DSYEV(), and vec.

Referenced by vec12().

                                                                {

    int i,j,k;
    char jobz, uplo;
    int n, lda, lwork, info;
    double *a, *w, *work;
#ifdef SR
        int *iwork, liwork;
        liwork = 5 * xNsize + 3;
        iwork = (int*)malloc(liwork * sizeof(double));
        lwork = 2 * xNsize * xNsize + 6 * xNsize + 1; /* 3*xNsize OK?*/
#else
        lwork = 4 * xNsize; /* 3*xNsize OK?*/
#endif

    n = lda = xNsize;

    a = (double*)malloc(xNsize*xNsize*sizeof(double));
    w = (double*)malloc(xNsize*sizeof(double));
    work = (double*)malloc(lwork*sizeof(double));

    k=0;
    for(j=0;j<xNsize;j++){
        for(i=0;i<xNsize;i++){
            a[k] = A[i][j];
            k++;
        }
    }

    jobz = 'V';
    uplo = 'U';

#ifdef SR
        M_DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, iwork, &liwork, &info);
        free(iwork);
#else
        dsyev_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
#endif

    k=0;
    for(i=0;i<xNsize;i++){
        for(j=0;j<xNsize;j++){
            vec[i][j]=a[k];
            k++;
        }
    }

    if(info != 0){
        free(a);
        free(w);
        free(work);
        return 0;
    }

    for(k=0;k<xNsize;k++){
        r[k] = w[k];
    }

    free(a);
    free(w);
    free(work);

    return 1;
}

DSEVXU()

int DSEVXU	(	int	xNsize,
		double **	A,
		double *	r,
		double **	X,
		int	xNev
	)

obtain eigenvalues A

Parameters

xNsize	size of A
A	matrix
r	eigenvalues
X	eigenvectors
xNev	number of eigenvalues

Returns

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 291 of file matrixlapack.c.

References dlamch_(), dsyevx_(), and X.

                                                                   {

    int i,j,k;
    char jobz, range, uplo;
    int n, lda, il, iu, m, ldz, lwork, *iwork, *ifail, info;
    double *a, *w, *work, vl, vu, abstol, *z;

    n = lda = ldz = xNsize;
    lwork = 8*xNsize;

    a = (double*)malloc(xNsize*xNsize*sizeof(double));
    w = (double*)malloc(xNsize*sizeof(double));
    z = (double*)malloc(xNsize*xNsize*sizeof(double));
    work  = (double*)malloc(lwork*sizeof(double));
    iwork = (int*)malloc(5*xNsize*sizeof(int));
    ifail = (int*)malloc(xNsize*sizeof(int));

#ifdef SR
        abstol = 0.0;
#else
        abstol = 2.0*dlamch_("S");
#endif
        vl = vu = 0.0;

    k=0;
    for(j=0;j<xNsize;j++){
        for(i=0;i<xNsize;i++){
            a[k] = A[i][j];
            k++;
        }
    }

    jobz  = 'V';
    range = 'I';
    uplo  = 'U';

    il = xNsize-xNev+1;
    iu = xNsize;
    m = iu-il+1;

    dsyevx_(&jobz, &range, &uplo, &n, a, &lda, &vl, &vu, &il, &iu, &abstol, 
            &m, w, z, &ldz, work, &lwork, iwork, ifail, &info);

    if(info != 0){
        free(a);
        free(w);
        free(z);
        free(work);
        free(iwork);
        free(ifail);
        return 0;
    }

    for(k=0;k<xNev;k++){
        r[k+xNsize-xNev] = w[k];
    }

    for(k=0;k<xNsize*xNev;k++){
        X[k%xNsize][k/xNsize+xNsize-xNev] = z[k];
    }
    free(a);
    free(w);
    free(z);
    free(work);
    free(iwork);
    free(ifail);

    return 1;
}

dsyev_()

int dsyev_	(	char *	jobz,
		char *	uplo,
		int *	n,
		double *	a,
		int *	lda,
		double *	w,
		double *	work,
		int *	lwork,
		int *	info
	)

Referenced by DSEVvector().

dsyevx_()

int dsyevx_	(	char *	jobz,
		char *	range,
		char *	uplo,
		int *	n,
		double *	a,
		int *	lda,
		double *	vl,
		double *	vu,
		int *	il,
		int *	iu,
		double *	abstol,
		int *	m,
		double *	w,
		double *	z__,
		int *	ldz,
		double *	work,
		int *	lwork,
		int *	iwork,
		int *	ifail,
		int *	info
	)

Referenced by DSEVXU().

M_DSYEV()

int M_DSYEV	(	char *	jobz,
		char *	uplo,
		int *	n,
		double *	a,
		int *	lda,
		double *	w,
		double *	work,
		int *	lwork,
		int *	info
	)

Referenced by DSEVvalue(), and DSEVvector().

to_f()

void to_f	(	int	N,
		int	M,
		double **	A,
		double *	a
	)

function for transforming Row-major matrix (C) to Column-major matrix (Fortran)

Parameters

[in]	N
[in]	M
[in]	A	Row-major matrix
[out]	a	Column-major matrix

Author

Takahiro Misawa (The University of Tokyo)

Parameters

[in]	N
[in]	M
[in]	A
[out]	a

Definition at line 67 of file matrixlapack.c.

Referenced by DSEVvalue().

 {
  int i,j,k;
  k=0;
  for(j=0;j<M;j++){
    for(i=0;i<N;i++){
      a[k] = A[i][j];
      k++;
    }
  }
}

zheev_()

int zheev_	(	char *	jobz,
		char *	uplo,
		int *	n,
		double complex *	a,
		int *	lda,
		double *	w,
		double complex *	work,
		int *	lwork,
		double *	rwork,
		int *	info
	)

Referenced by ZHEEVall().

ZHEEVall()

int ZHEEVall	(	int	xNsize,
		double complex **	A,
		double complex *	r,
		double complex **	vec
	)

obtain eigenvalues and eigenvectors of Hermite matrix A

Parameters

xNsize	size of matrix
A	matrix
r	eigenvalues
vec	eigenvectors

Returns

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 375 of file matrixlapack.c.

References vec, zheev_(), and zheevd_().

Referenced by lapack_diag().

                                                                                    {

    int i,j,k;
    char jobz, uplo;
    int n, lda, lwork, info, lrwork;
        double *rwork;
    double *w;
    double complex *a, *work;
#ifdef SR
        int *iwork, liwork;
        liwork = 5 * xNsize + 3;
        iwork = (int*)malloc(liwork * sizeof(double));
#endif

    n = lda = xNsize;
#ifdef SR
        lwork = xNsize*xNsize + 2 * xNsize; /* 3*xNsize OK?*/
        lrwork = 2 * xNsize*xNsize + 5*xNsize + 1;
#else
        lwork = 4*xNsize; /* 3*xNsize OK?*/
        lrwork = lwork;
#endif

    a = (double complex*)malloc(xNsize*xNsize*sizeof(double complex));
    w = (double*)malloc(xNsize*sizeof(double));
    work = (double complex*)malloc(lwork*sizeof(double complex));
    rwork = (double*)malloc(lrwork*sizeof(double));

    k=0;
    for(j=0;j<xNsize;j++){
        for(i=0;i<xNsize;i++){
            a[k] = A[i][j];
            k++;
        }
    }

    jobz = 'V';
    uplo = 'U';

#ifdef SR
        int zheevd_(char *jobz, char *uplo, int *n, double complex *a, int *lda, double *w, double complex *work, int *lwork, double *rwork, int *iwork, int *liwork, int *info); 
        free(iwork);
#else
    zheev_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, rwork, &info);
#endif

    if(info != 0){
        free(a);
        free(w);
        free(work);
        free(rwork);
        return 0;
    }

    k=0;
    for(i=0;i<xNsize;i++){
        for(j=0;j<xNsize;j++){
            vec[i][j]=a[k];
            k++;
        }
    }

  
    for(k=0;k<xNsize;k++){
        r[k] = w[k];
    }

    free(a);
    free(w);
    free(work);
    free(rwork);

    return 1;
}