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 "common/setmemory.h"
#include "xsetmem.h"

Include dependency graph for vec12.c:

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 BindStruct::Def, DSEVvector(), DefineList::iCalcType, DefineList::iFlgCalcSpec, DefineList::LanczosTarget, DefineList::nvec, stdoutMPI, and vec.

Referenced by CalcSpectrumByTPQ(), and Lanczos_EigenValue().

  {
  unsigned int j,k,nvec;

  double **tmpA, **tmpvec;
  double *tmpr;
    
  nvec = X->Def.nvec;
  tmpA = d_2d_allocate(ndim, ndim);
  tmpvec = d_2d_allocate(ndim, ndim);
  tmpr = d_1d_allocate(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_d_2d_allocate(tmpA);
  free_d_1d_allocate(tmpr);
  free_d_2d_allocate(tmpvec);
}/*void vec12*/

Here is the call graph for this function:

Here is the caller graph for this function: