vec12.c File Reference

Functions to Diagonalize a tri-diagonal matrix and store eigenvectors into vec. More...

#include "matrixlapack.h"
#include "Common.h"
#include "wrapperMPI.h"
#include "mfmemory.h"
#include "xsetmem.h"

Go to the source code of this file.

Functions
void	vec12 (double alpha[], double beta[], unsigned int ndim, double tmp_E[], struct BindStruct *X)
	Diagonalize a tri-diagonal matrix and store eigenvectors into vec. More...

Detailed Description

Functions to Diagonalize a tri-diagonal matrix and store eigenvectors into vec.

Definition in file vec12.c.

Function Documentation

vec12()

void vec12	(	double	alpha[],
		double	beta[],
		unsigned int	ndim,
		double	tmp_E[],
		struct BindStruct *	X
	)

Diagonalize a tri-diagonal matrix and store eigenvectors into vec.

Author

Takahiro Misawa (The University of Tokyo)

Youhei Yamaji (The University of Tokyo)

Definition at line 31 of file vec12.c.

References alpha, beta, DSEVvector(), stdoutMPI, vec, and X.

Referenced by CalcSpectrumByTPQ(), and Lanczos_EigenValue().

  {
  unsigned int j,k,nvec;

  double **tmpA, **tmpvec;
  double *tmpr;
    
  nvec = X->Def.nvec;
  d_malloc2(tmpA,ndim,ndim); 
  d_malloc2(tmpvec,ndim,ndim); 
  d_malloc1(tmpr,ndim);

#pragma omp parallel for default(none) firstprivate(ndim) private(j,k) shared(tmpA)
  for(k=0;k<=ndim-1;k++)
    for(j=0;j<=ndim-1;j++) tmpA[k][j]=0.0;
#pragma omp parallel for default(none) firstprivate(ndim, nvec) private(j,k) shared(vec)
  for(k=1;k<=nvec;k++)
    for(j=1;j<=ndim;j++) vec[k][j]=0.0;

#pragma omp parallel for default(none) firstprivate(ndim, alpha, beta) private(j) shared(tmpA)
  for(j=0;j<=ndim-2;j++){
    tmpA[j][j]=alpha[j+1];
    tmpA[j][j+1]=beta[j+1];
    tmpA[j+1][j]=beta[j+1];
  }/*for(j=0;j<=ndim-2;j++)*/
  tmpA[ndim-1][ndim-1]=alpha[ndim];

  DSEVvector( ndim, tmpA, tmpr, tmpvec );
  if(X->Def.iCalcType==Lanczos && X->Def.iFlgCalcSpec == 0)
    fprintf(stdoutMPI, "  Lanczos EigenValue in vec12 = %.10lf \n ",tmpr[0]);
 
  if (nvec <= ndim) {
    if (nvec < X->Def.LanczosTarget) nvec = X->Def.LanczosTarget;
    
#pragma omp parallel for default(none) firstprivate(ndim, nvec) private(j,k) shared(tmpvec, vec, tmp_E, tmpr)
    for(k=1;k<=nvec;k++){
      tmp_E[k]=tmpr[k-1];
      for (j = 1; j <= ndim; j++) vec[k][j] = tmpvec[k - 1][j - 1];
    }/*for(k=1;k<=nvec;k++)*/
  }/*if(nvec<=ndim)*/
  else{
#pragma omp parallel for default(none) firstprivate(ndim, nvec) private(j,k) shared(tmpvec, vec, tmp_E, tmpr)
    for(k=1;k<=ndim;k++){
      tmp_E[k]=tmpr[k-1];
      for (j = 1; j <= ndim; j++) vec[k][j] = tmpvec[k - 1][j - 1];
    }/*for(k=1;k<=ndim;k++)*/
  }/*if(nvec>ndim)*/
  free(tmpA);
  free(tmpr);
  free(tmpvec);
}/*void vec12*/