matrixlapack.c File Reference

wrapper for linear algebra operations using lapack More...

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

Include dependency graph for matrixlapack.c:

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().

Here is the caller graph for this function:

dgetri_()

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

Referenced by DInv().

Here is the caller graph for this function:

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 153 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;
}

Here is the call graph for this function:

dlamch_()

double dlamch_

(

char *

cmach

)

Referenced by DSEVXU().

Here is the caller graph for this function:

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 93 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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 209 of file matrixlapack.c.

References dsyev_(), and M_DSYEV().

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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 290 of file matrixlapack.c.

References dlamch_(), and dsyevx_().

                                                                   {

  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;
}

Here is the call graph for this function:

dsyev_()

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

Referenced by DSEVvector().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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 66 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++;
    }
  }
}

Here is the caller graph for this function:

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().

Here is the caller graph for this function:

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 374 of file matrixlapack.c.

References 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;
}

Here is the call graph for this function:

Here is the caller graph for this function: