matrixlapack.cpp File Reference

wrapper for linear algebra operations using lapack More...

#include "matrixlapack.hpp"
#include <cstdlib>

Go to the source code of this file.

Functions
int	zheev_ (char *jobz, char *uplo, int *n, std::complex< double > *a, int *lda, double *w, std::complex< double > *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)
int	ZHEEVall (int xNsize, std::complex< double > **A, double *r, std::complex< double > **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.cpp.

Function Documentation

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
	)

zheev_()

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

Referenced by ZHEEVall().

ZHEEVall()

int ZHEEVall	(	int	xNsize,
		std::complex< double > **	A,
		double *	r,
		std::complex< double > **	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 63 of file matrixlapack.cpp.

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;
  std::complex<double> *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 = (std::complex<double>*)malloc(xNsize*xNsize * sizeof(std::complex<double>));
  w = (double*)malloc(xNsize * sizeof(double));
  work = (std::complex<double>*)malloc(lwork * sizeof(std::complex<double>));
  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, std::complex<double> *a, int *lda, double *w, std::complex<double> *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[j][i] = a[k];
      k++;
    }
  }


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

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

  return 1;
}