IPCC  1.0
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
KNLanczosMethod.cpp
Go to the documentation of this file.
1 
7 #include "stdafx.h"
8 #include "KNLanczosMethod.h"
9 #include "KNTimeMeasurement.h"
10 #include <math.h>
11 #include <string.h>
12 #include <algorithm>
13 #include "CKNGlobal.h"
14 #include "KNIPCCUtility.h"
15 #include "KNMPIManager.h"
16 #include "mkl.h"
17 
18 #ifdef _WIN32
19 #include <direct.h>
20 #else
21 #include <sys/stat.h>
22 #include <sys/types.h>
23 #endif
24 
25 
26 #include "XeonPhi_header.h"
27 unsigned int X_largest;
28 
29 using namespace std;
30 bool CKNLanczosMethod::m_bStop = false;
31 
32 CKNLanczosMethod::CKNLanczosMethod()
33 {
34  m_pV = NULL;
35  m_pAMatrix = NULL;
36  m_pAMyLocalBlock = NULL;
37  m_pALeftBlock = NULL;
38  m_pARightBoloc = NULL;
39  m_nMatrixSize = 0;
40  m_nIterationCount = 0;
41  m_nEigenValueCheckInterval = 0;
42  m_fEigenvalueMin = 0;
43  m_fEignevalueMax = 0;
44  m_bReorthogonalization = false;
45  m_bCalcuEigenvector = false;
46  //m_nPrevIterationCount = 0;
47 
48  m_pEigenValues = NULL;
49  m_pEigenVectors = NULL;
50  m_pConvergedEigenValues = NULL;
51  m_pConvergedEigenVectors = NULL;
52  m_pRangeCheckedEigenValues = NULL;
53  m_pRangeCheckedEigenVectors = NULL;
54  m_pNoneSpuriousValues = NULL;
55  m_pNoneSpuriousVectors = NULL;
56  m_pNonClustersValues = NULL;
57  m_pNonClustersVectors = NULL;
58 
59  m_pRangecheckedIndex = NULL;
60  m_pNonSpuriousValueIndex = NULL;
61  m_pConvergedIndex = NULL;
62  m_pNonClustersValueIndex = NULL;
63  m_pCheckNonClusterValue = NULL;
64  m_floadMIC = 0.0;
65 }
66 
67 CKNLanczosMethod::~CKNLanczosMethod()
68 {
69  FinalizeTemporaryArrayAndVector();
70 }
71 
76 bool CKNLanczosMethod::InitializeTemporaryArrayAndVector(int nIterationCount)
77 {
78  bool bRtn = true;
79  /*if (nIterationCount == m_nPrevIterationCount)
80  return bRtn;*/
81 
82  //if (0 != m_nPrevIterationCount)
83  FinalizeTemporaryArrayAndVector();
84 
85  bRtn = false;
86 
87  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
88  ALLOC_WITH_NULL_INIT(m_pEigenValues, double, nIterationCount);
89  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
90 
91  bRtn = true;
92  return bRtn;
93 }
94 
95 void CKNLanczosMethod::FinalizeTemporaryArrayAndVector()
96 {
97  FREE_MEM(m_pEigenValues);
98  FREE_MEM(m_pEigenVectors);
99  FREE_MEM(m_pConvergedEigenValues);
100  FREE_MEM(m_pConvergedEigenVectors);
101  FREE_MEM(m_pRangeCheckedEigenValues);
102  FREE_MEM(m_pRangeCheckedEigenVectors);
103  FREE_MEM(m_pNoneSpuriousValues);
104  FREE_MEM(m_pNoneSpuriousVectors);
105  FREE_MEM(m_pNonClustersValues);
106  FREE_MEM(m_pNonClustersVectors);
107  FREE_MEM(m_pRangecheckedIndex);
108  FREE_MEM(m_pNonSpuriousValueIndex);
109  FREE_MEM(m_pConvergedIndex);
110  FREE_MEM(m_pNonClustersValueIndex);
111  FREE_MEM(m_pCheckNonClusterValue);
112 }
113 
130 CKNLanczosMethod::LPEIGENVALUE_RESULT CKNLanczosMethod::DoLanczosMethod(CKNMatrixOperation::CKNCSR *pAMatrix, unsigned int nIterationCount, unsigned int nEigenValueCheckInterval, unsigned int nEigenValueCount, double fEigenvalueMin, double fEignevalueMax, double fConvergenceTolerance, bool bReorthogonalization, bool bCalcuEigVector, bool bWaveFunction, double load_in_MIC, CKNMatrixOperation::CKNCSR *pmylocalblock, CKNMatrixOperation::CKNCSR *leftlocalblock, CKNMatrixOperation::CKNCSR *rightlocalblock)
131 {
132  LPEIGENVALUE_RESULT lpRtn = NULL;
133 
134  if (NULL == pAMatrix)
135  {
136  throw ERROR_MALLOC;
137  return NULL;
138  }
139 
140  InitVariables();
141 
142  m_pAMatrix = pAMatrix;
143  m_pAMyLocalBlock = pmylocalblock;
144  m_pALeftBlock = leftlocalblock;
145  m_pARightBoloc = rightlocalblock;
146  m_nMatrixSize = pAMatrix->GetColumnCount();
147  m_nIterationCount = nIterationCount;
148  m_nEigenValueCheckInterval = nEigenValueCheckInterval;
149  m_nEigenValueCount = nEigenValueCount;
150  m_fEigenvalueMin = fEigenvalueMin;
151  m_fEignevalueMax = fEignevalueMax;
152  m_bReorthogonalization = bReorthogonalization;
153  m_bCalcuEigenvector = bCalcuEigVector;
154  m_fConvergenceTolerance = fConvergenceTolerance;
155  m_floadMIC = load_in_MIC;
156 
157  lpRtn = LanczosIteration();
158  if( bCalcuEigVector )
159  DoResidualCheck(pAMatrix, lpRtn);
160 
161  //if (bWaveFunction && CKNMPIManager::IsRootRank())
162  if (bWaveFunction)
163  BuildWaveFunction(lpRtn);
164 
165  //m_nPrevIterationCount = nIterationCount;
166 
167  return lpRtn;
168 }
169 
173 CKNLanczosMethod::LPEIGENVALUE_RESULT CKNLanczosMethod::LanczosIteration()
174 {
175  CKNComplex *pAlpha = NULL;
176  double *pAlphaReal = NULL;
177  double *pBeta = NULL;
178  double *pWj = NULL;
179  double *pWjm1 = NULL;
180  double *pWjp1 = NULL;
181  CKNMatrixOperation::CKNVector *pW = NULL;
182  CKNMatrixOperation::CKNVector V1;
183  double *pCalculatedEigenVector = NULL;
184  unsigned int i;
185 
186  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
187  LPEIGENVALUE_RESULT lpResult = (LPEIGENVALUE_RESULT)malloc(sizeof(EIGENVALUE_RESULT));
188  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
189 
190  InitLanczosIterationVariables(&pAlpha, &pAlphaReal, &pBeta, &pWj, &pWjm1, &pWjp1, &pW);
191  InitLanczosVector();
192 
193  pAlphaReal[0] = pBeta[0] = pBeta[1] = 0;
194 
195 #ifdef DISABLE_MPI_ROUTINE
196  V1.SetSize(m_nMatrixSize);
197  V1.SetAt(0, 1, 0);
198 #else //DISABLE_MPI_ROUTINE
199  V1.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
200  if( CKNMPIManager::IsMultiLevelMPI())
201  {
202  /*V1.BuildRandomVector();
203  V1.Normalize(true);*/
204  if (CKNMPIManager::IsRootRank())
205  V1.SetAt(CKNMPIManager::GetLanczosGroupIndex()*2, 1, 0);
206  }
207  else
208  {
209  if (CKNMPIManager::IsRootRank())
210  V1.SetAt(0, 1, 0);
211  }
212 #endif //DISABLE_MPI_ROUTINE
213 
214  lpResult->nEigenValueCount = 0;
215  lpResult->pEigenValues = NULL;
216  lpResult->pDegeneratedIndex = NULL;
217  lpResult->pEigenVectors = NULL;
218  lpResult->nEigenValueCount = 0;
219  lpResult->pEigenValueFoundIteration = NULL;
220  lpResult->pEigenVectorsForAMatrix = NULL;
221  lpResult->pWaveFunctions = NULL;
222  lpResult->nMaxEigenValueFoundIteration = 0;
223  lpResult->nDegeneratedEigenValueCount = 0;
224 
225  if (m_bReorthogonalization && CKNMPIManager::IsRootRank())
226  {
227  m_pV[1] = V1;
228 
229  pWj[2] = 0;
230  pWj[3] = 0;
231  }
232 
234  LanczosIterationLoop(lpResult, &V1, m_nIterationCount, pAlpha, pAlphaReal, pBeta, pWj, pWjm1, pWjp1);
235  if (CKNLanczosMethod::m_bStop)
236  {
237  V1.Finalize();
238  FinalLanczosVector();
239  FinalizeLanczosInterationVariable(pAlpha, pAlphaReal, pBeta, pWj, pWjm1, pWjp1, pW);
240 
241  return lpResult;
242  }
243 
245  if (m_bCalcuEigenvector && false == CKNLanczosMethod::m_bStop)
246  {
247 
248 #ifndef DISABLE_MPI_ROUTINE
249  CKNMPIManager::BroadcastLanczosResult(lpResult, lpResult->nMaxEigenValueFoundIteration);
250 #endif
251  if (lpResult->nEigenValueCount > 0)
252  {
253  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
254  lpResult->pEigenVectorsForAMatrix = new CKNMatrixOperation::CKNVector[lpResult->nEigenValueCount];
255  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
256  for (i = 0; i < lpResult->nEigenValueCount; i++)
257 #ifdef DISABLE_MPI_ROUTINE
258  lpResult->pEigenVectorsForAMatrix[i].SetSize(m_nMatrixSize);
259 #else //DISABLE_MPI_ROUTINE
260  lpResult->pEigenVectorsForAMatrix[i].SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
261 #endif
262 
263  if (m_bReorthogonalization)
264  {
265  for (i = 1; i <= lpResult->nMaxEigenValueFoundIteration; i++)
266  CalculateEigenVector(lpResult, m_pV[i], i);
267  }
268  else
269  LanczosIterationLoop(lpResult, &V1, lpResult->nMaxEigenValueFoundIteration, pAlpha, pAlphaReal, pBeta, pWj, pWjm1, pWjp1, true);
270  }
271 
272  if (CKNLanczosMethod::m_bStop)
273  return lpResult;
274 
275  for (i = 0; i < lpResult->nEigenValueCount; ++i)
276 #ifndef DISABLE_MPI_ROUTINE
277  lpResult->pEigenVectorsForAMatrix[i].Normalize(true);
278 #else //DISABLE_MPI_ROUTINE
279  lpResult->pEigenVectorsForAMatrix[i].Normalize();
280 #endif //DISABLE_MPI_ROUTINE
281  }
282 
283  V1.Finalize();
284 
285  FinalLanczosVector();
286  FinalizeLanczosInterationVariable(pAlpha, pAlphaReal, pBeta, pWj, pWjm1, pWjp1, pW);
287 
288  return lpResult;
289 }
290 
303 void CKNLanczosMethod::LanczosIterationLoop(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, CKNMatrixOperation::CKNVector *V1, unsigned int nIterationCount, CKNComplex *pAlpha, double *pAlphaReal, double *pBeta, double *pWj, double *pWjm1, double *pWjp1, bool bMakeEigvVector)
304 {
305 // CKNMatrixOperation::CKNVector vector1, vector2;
306  CKNMatrixOperation::CKNVector Vj, Vjp1, Vjm1, W, VTemp;;
307  double fANorm = 0;
308  int nEigenvalueSolvingCount = 0;
309  int nEigenvalueSolvingFinal = m_nIterationCount / m_nEigenValueCheckInterval;
310  unsigned int j;
311  unsigned int nSizePHI;
312  char szMsg[1024];
313 
314 
315  Vj = *V1;
316 
317 #ifdef DISABLE_MPI_ROUTINE
318  Vjp1 = *V1;
319  Vjm1 = *V1;
320  W = *V1;
321 #else
322  Vjp1.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
323  Vjm1.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
324  W.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
325 
326 
327  int nPrevRank = (CKNMPIManager::GetCurrentRank() - 1 + CKNMPIManager::GetTotalNodeCount()) % CKNMPIManager::GetTotalNodeCount();
328  int nNextRank = (CKNMPIManager::GetCurrentRank() + 1) % CKNMPIManager::GetTotalNodeCount();
329  int tag = 1002;
330  unsigned int nSizeFromPrevRank, nSizeFromNextRank;
331  unsigned int nSizetoPrevRank, nSizetoNextRank;
332  MPI_Status stat_sr[2];
333  MPI_Request req_sr[2];
334 
335  nSizetoPrevRank = m_pAMatrix->nComponentsFirstUnitCell;
336  nSizetoNextRank = m_pAMatrix->nComponentsLastUnitCell;
337 
338  //printf("Rank %d: Preparing to get sizes. PrevRank %d, NextRank %d\n", CKNMPIManager::GetCurrentRank(), nPrevRank, nNextRank);
339 
340  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::COMM);
341  MPI_Irecv(&nSizeFromPrevRank, 1, MPI_INT, nPrevRank, tag, CKNMPIManager::GetMPIComm(), &req_sr[0]);
342  MPI_Isend(&nSizetoNextRank, 1, MPI_INT, nNextRank, tag, CKNMPIManager::GetMPIComm(), &req_sr[1]);
343  MPI_Waitall(2, req_sr, stat_sr); // now Brt has B of right neighbor
344 
345  MPI_Irecv(&nSizeFromNextRank, 1, MPI_INT, nNextRank, tag, CKNMPIManager::GetMPIComm(), &req_sr[0]);
346  MPI_Isend(&nSizetoPrevRank, 1, MPI_INT, nPrevRank, tag, CKNMPIManager::GetMPIComm(), &req_sr[1]);
347  MPI_Waitall(2, req_sr, stat_sr);
348  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::COMM);
349 
350  //printf("Rank %d: P=%d, M=%d, N=%d\n", CKNMPIManager::GetCurrentRank(), nSizeFromPrevRank, CKNMPIManager::GetCurrentLoadBalanceCount(), nSizeFromNextRank);
351 
352 
353 #endif
354 
355  nSizePHI = (int)(((double)CKNMPIManager::GetCurrentLoadBalanceCount()) * m_floadMIC / 100.0);
356  m_floadMIC = ((double)nSizePHI) / ((double)CKNMPIManager::GetCurrentLoadBalanceCount())*100.0;
357 
358  sprintf(szMsg, "-[Rank %03d] MIC-load adjusted to %.1f(%%). DOF(MIC, Total) = (%d, %d)\n", CKNMPIManager::GetCurrentRank(), m_floadMIC, nSizePHI, CKNMPIManager::GetCurrentLoadBalanceCount());
359  CKNIPCCUtility::ShowMsg(szMsg);
360 
361  double *input_real = Vj.m_vectValueRealBuffer.data();
362  double *input_imaginary = Vj.m_vectValueImaginaryBuffer.data();
363  unsigned int input_real_size = Vj.m_vectValueRealBuffer.size();
364  unsigned int input_imaginary_size = Vj.m_vectValueImaginaryBuffer.size();
365  //#pragma offload_transfer target(mic:phi_tid) nocopy(input_real[0:input_real_size] : ALLOC)
366  //#pragma offload_transfer target(mic:phi_tid) nocopy(input_imaginary[0:input_imaginary_size] : ALLOC)
367 
368  double *output_real = W.m_vectValueRealBuffer.data();
369  double *output_imaginary = W.m_vectValueImaginaryBuffer.data();
370  //unsigned int output_real_size = W.m_vectValueRealBuffer.size();
371  //unsigned int output_imaginary_size = W.m_vectValueImaginaryBuffer.size();
372  unsigned int output_real_size = nSizePHI;
373  unsigned int output_imaginary_size = nSizePHI;
374 
375  if( 0 == nSizePHI )
376  {
377  output_real_size = 1;
378  output_imaginary_size = 1;
379  }
380 
381 #pragma offload_transfer target(mic:phi_tid) nocopy(output_real[0:output_real_size] : ALLOC)
382 #pragma offload_transfer target(mic:phi_tid) nocopy(output_imaginary[0:output_imaginary_size] : ALLOC)
383  // if(CKNMPIManager::GetTotalNodeCount() <= 3)
384  VTemp.SetSize(m_nMatrixSize);
385  // else
386  //VTemp.SetSize(nSizeFromPrevRank+nSizeFromNextRank+CKNMPIManager::GetCurrentLoadBalanceCount());
387 
388  double *vtemp_real = VTemp.m_vectValueRealBuffer.data();
389  double *vtemp_imaginary = VTemp.m_vectValueImaginaryBuffer.data();
390  unsigned int vtemp_size = VTemp.m_vectValueRealBuffer.size();
391 
392 #pragma offload_transfer target(mic:phi_tid) nocopy(vtemp_real[0:vtemp_size] : ALLOC)
393 #pragma offload_transfer target(mic:phi_tid) nocopy(vtemp_imaginary[0:vtemp_size] : ALLOC)
394 #pragma offload_transfer target(mic:phi_tid) in(vtemp_real[0:vtemp_size] : REUSE)
395 #pragma offload_transfer target(mic:phi_tid) in(vtemp_imaginary[0:vtemp_size] : REUSE)
396  // X_largest = CKNMPIManager::GetLoadBalanceCount(0);
397  // for(int i = 1; i < CKNMPIManager::GetTotalNodeCount(); i++)
398  // {
399  // unsigned int Btemp = CKNMPIManager::GetLoadBalanceCount(i);
400  // if(Btemp > X_largest)
401  // X_largest = Btemp;
402  // }
403 
404  // CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
405  // double *X = (double *)_mm_malloc(X_largest * 2 * sizeof(double), 64);
406  // double *Xrt = (double *)_mm_malloc(X_largest * 2 * sizeof(double), 64);
407  // double *Xlt = (double *)_mm_malloc(X_largest * 2 * sizeof(double), 64);
408  //#pragma offload_transfer target(mic:phi_tid) nocopy(X[0:X_largest * 2] : align(64) ALLOC)
409  //#pragma offload_transfer target(mic:phi_tid) nocopy(Xrt[0:X_largest * 2] : align(64) ALLOC)
410  //#pragma offload_transfer target(mic:phi_tid) nocopy(Xlt[0:X_largest * 2] : align(64) ALLOC)
411  // double *X = new double[X_largest*2];
412  // double *Xrt = new double[X_largest*2];
413  // double *Xlt = new double[X_largest*2];
414  //#pragma offload_transfer target(mic:phi_tid) nocopy(X[0:X_largest * 2] : ALLOC)
415  //#pragma offload_transfer target(mic:phi_tid) nocopy(Xrt[0:X_largest * 2] : ALLOC)
416  //#pragma offload_transfer target(mic:phi_tid) nocopy(Xlt[0:X_largest * 2] : ALLOC)
417  // CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
418 
419  char szBuffer[1024];
420  int nSize;
421 #ifdef DISABLE_MPI_ROUTINE
422  nSize = m_nMatrixSize;
423 #else
424  nSize = CKNMPIManager::GetCurrentLoadBalanceCount();
425 #endif
426 
427  for (j = 1; j <= nIterationCount; j++)
428  {
429 
430  if (0 == j % 500 && CKNMPIManager::IsRootRank())
431  {
432  sprintf(szMsg, "[#%8d] Lanczos interation going on\n", j);
433  CKNIPCCUtility::ShowMsg(szMsg);
434  }
435 
436  if (bMakeEigvVector)
437  CalculateEigenVector(lpResult, Vj, j);
438 
439 #ifdef DISABLE_MPI_ROUTINE
440  if (j == 1)
441  CKNIPCCUtility::ShowMsg("-Using MVMulOptimal with no offload\n");
442  CKNMatrixOperation::MVMulOptimal(m_pAMatrix, &Vj, &W);
443 #else //DISABLE_MPI_ROUTINE
444  /* /// Using Normal MV Mul
445  if (j == 1)
446  if(CKNMPIManager::IsRootRank())
447  CKNIPCCUtility::ShowMsg("-Using MVMul\n");
448  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MVMUL);
449  CKNMatrixOperation::MVMul(m_pAMatrix, &Vj, &W);
450  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MVMUL);
451  */
452 
453  /* /// Using Async block communication
454  if(CKNMPIManager::GetTotalNodeCount() <= 2)
455  {
456  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MVMUL);
457  CKNMatrixOperation::MVMulOptimal(m_pAMatrix, &Vj, &W); /// wj <- Avj
458  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MVMUL);
459  if (j == 1)
460  if(CKNMPIManager::IsRootRank())
461  CKNIPCCUtility::ShowMsg("-Using MVMulOptimal\n");
462  }
463  else
464  {
465  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MVMUL);
466  CKNMatrixOperation::MVMulEx_AsyncCommWithLocalBlocks(m_pAMyLocalBlock, m_pALeftBlock, m_pARightBoloc, &Vj, &W, X, Xrt, Xlt);
467  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MVMUL);
468  if (j == 1)
469  if(CKNMPIManager::IsRootRank())
470  CKNIPCCUtility::ShowMsg("-Using MVMulEx_AsyncCommWithLocalBlocks\n");
471  }
472  */
474  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MVMUL);
475  CKNMatrixOperation::MVMulEx_Optimal(m_pAMatrix, &Vj, &W, nSizeFromPrevRank, nSizeFromNextRank, &VTemp, nSizePHI);
476  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MVMUL);
477  if (j == 1)
478  if (CKNMPIManager::IsRootRank())
479  printf("-Using MVMulEx_Optimal with %.1f(%%) offload\n", m_floadMIC);
480 #endif //DISABLE_MPI_ROUTINE
481 
482  if (1 != j)
483  W.ScalarMultiThanMinusVector(pBeta[j], &Vjm1);
484 
485  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::VVDOT);
486  CKNMatrixOperation::VVDot(&W, &Vj, &pAlpha[j]);
487  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::VVDOT);
488  pAlphaReal[j] = pAlpha[j].GetRealNumber();
489  W.ScalarMultiThanMinusVector(pAlphaReal[j], &Vj);
490 
491 #ifdef DISABLE_MPI_ROUTINE
492  pBeta[j + 1] = W.GetNorm();
493 #else
494  pBeta[j + 1] = W.GetNorm(true);
495 #endif
496  Vjp1 = W;
497  Vjp1.ScalarDivision(pBeta[j + 1]);
498 
499  if (m_bReorthogonalization && false == bMakeEigvVector)
500  {
501  if (CKNMPIManager::IsRootRank())
502  {
503  if (1 == j && fANorm < fabs(pAlphaReal[1]) + pBeta[2])
504  fANorm = fabs(pAlphaReal[1]) + pBeta[2];
505  else if (fANorm < fabs(pAlphaReal[j]) + pBeta[j + 1] + pBeta[j])
506  fANorm = fabs(pAlphaReal[j]) + pBeta[j + 1] + pBeta[j];
507  }
508 
509  m_pV[j + 1] = Vjp1;
510  if (CheckAndDoSelectiveReorthogonalization(j-1, pAlphaReal+1, pBeta+2, pWj, pWjm1, pWjp1, fANorm))
511  {
512  Vj = m_pV[j];
513  Vjp1 = m_pV[j + 1];
514  }
515  }
516 
517  if (CKNLanczosMethod::m_bStop)
518  break;
519 
520  if (0 == j % m_nEigenValueCheckInterval && false == bMakeEigvVector)
521  {
522  bool bFound = false;
523  int nPrevEVCount = lpResult->nEigenValueCount;
524 
525  if (InitializeTemporaryArrayAndVector(j))
526  {
527  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::EVALUE);
528  bFound = DoEigenValueSolving(j, pAlphaReal, pBeta, fANorm, lpResult, ++nEigenvalueSolvingCount == nEigenvalueSolvingFinal);
529  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::EVALUE);
530  }
531  else
532  {
533  bFound = true; // Can't alloc memeory, temporary code
534  CKNIPCCUtility::ShowMsg("\n[Warning] For memory allocation bug, can't calculate eigenvalue anymore!.\n\n");
535  }
536 
537  FinalizeTemporaryArrayAndVector();
538 
539  CKNMPIManager::BroadcastBool(&bFound);
540  if (bFound)
541  break;
542  }
543 
544  Vjm1 = Vj;
545  Vj = Vjp1;
546  }
547 
548  //#pragma offload_transfer target(mic:phi_tid) nocopy(input_real : FREE)
549  //#pragma offload_transfer target(mic:phi_tid) nocopy(input_imaginary : FREE)
550 
551 #pragma offload_transfer target(mic:phi_tid) nocopy(output_real : FREE)
552 #pragma offload_transfer target(mic:phi_tid) nocopy(output_imaginary : FREE)
553 
554 #pragma offload_transfer target(mic:phi_tid) nocopy(vtemp_real : FREE)
555 #pragma offload_transfer target(mic:phi_tid) nocopy(vtemp_imaginary : FREE)
556 
557  // CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::FREE_MEM);
558  // _mm_free(X);
559  // _mm_free(Xrt);
560  // _mm_free(Xlt);
561  // delete [] X;
562  // delete [] Xrt;
563  // delete [] Xlt;
564  //#pragma offload_transfer target(mic:phi_tid) nocopy(X : FREE)
565  //#pragma offload_transfer target(mic:phi_tid) nocopy(Xrt : FREE)
566  //#pragma offload_transfer target(mic:phi_tid) nocopy(Xlt : FREE)
567  // CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::FREE_MEM);
568 
569  Vj.Finalize();
570  Vjp1.Finalize();
571  Vjm1.Finalize();
572  W.Finalize();
573 }
574 
575 
576 void CKNLanczosMethod::InitLanczosVector()
577 {
578  unsigned int i;
579 
580  if (false == m_bReorthogonalization)
581  return;
582 
583  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
584  m_pV = new CKNMatrixOperation::CKNVector[m_nIterationCount + 2];
585  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
586  for (i = 0; i < m_nIterationCount + 2; i++)
587 #ifdef DISABLE_MPI_ROUTINE
588  m_pV[i].SetSize(m_nMatrixSize);
589 #else //DISABLE_MPI_ROUTINE
590  m_pV[i].SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
591 #endif //DISABLE_MPI_ROUTINE
592 
593 }
594 
595 void CKNLanczosMethod::FinalLanczosVector()
596 {
597  unsigned int i;
598  if (false == m_bReorthogonalization)
599  return;
600  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::FREE_MEM);
601  for (i = 0; i < m_nIterationCount + 2; i++)
602  m_pV[i].Finalize();
603 
604  delete[] m_pV;
605  m_pV = NULL;
606  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::FREE_MEM);
607 }
608 
617 void CKNLanczosMethod::InitLanczosIterationVariables(CKNComplex **pAlpha, double **pAlphaReal, double **pBeta, double **pWj, double **pWjm1, double **pWjp1, CKNMatrixOperation::CKNVector **pW)
618 {
619  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
620  *pAlpha = new CKNComplex[m_nIterationCount + 2];
621  *pAlphaReal = (double*)malloc(sizeof(double)*(m_nIterationCount + 5));
622  *pBeta = (double*)malloc(sizeof(double)*(m_nIterationCount + 5));
623 
624  if (m_bReorthogonalization && CKNMPIManager::IsRootRank())
625  {
626  *pWj = (double*)malloc(sizeof(double)*(m_nIterationCount + 3));
627  *pWjm1 = (double*)malloc(sizeof(double)*(m_nIterationCount + 3));
628  *pWjp1 = (double*)malloc(sizeof(double)*(m_nIterationCount + 3));
629 
630  memset(*pWj, 0, sizeof(double)*(m_nIterationCount + 3));
631  memset(*pWjm1, 0, sizeof(double)*(m_nIterationCount + 3));
632  memset(*pWjp1, 0, sizeof(double)*(m_nIterationCount + 3));
633  }
634  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
635 }
636 
645 void CKNLanczosMethod::FinalizeLanczosInterationVariable(CKNComplex *pAlpha, double *pAlphaReal, double *pBeta, double *pWj, double *pWjm1, double *pWjp1, CKNMatrixOperation::CKNVector *pW)
646 {
647  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::FREE_MEM);
648  if (m_bReorthogonalization && CKNMPIManager::IsRootRank())
649  {
650  FREE_MEM(pWj);
651  FREE_MEM(pWjm1);
652  FREE_MEM(pWjp1);
653  }
654 
655  delete[] pAlpha;
656  pAlpha = NULL;
657 
658  free(pAlphaReal);
659  pAlphaReal = NULL;
660 
661  free(pBeta);
662  pBeta = NULL;
663  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::FREE_MEM);
664 }
665 
666 
667 
673 void CKNLanczosMethod::CalculateEigenVector(LPEIGENVALUE_RESULT lpResult, CKNMatrixOperation::CKNVector V, unsigned int nIterationIndex)
674 {
675  unsigned int i, k, nRepeatCount;
676 
678  for (i = 0; i < lpResult->nEigenValueCount; i++)
679  {
680  if (nIterationIndex > lpResult->pEigenValueFoundIteration[i])
681  continue;
682 
683 #ifdef DISABLE_MPI_ROUTINE
684  nRepeatCount = m_nMatrixSize;
685 #else //DISABLE_MPI_ROUTINE
686  nRepeatCount = CKNMPIManager::GetCurrentLoadBalanceCount();
687 #endif
688 
689  for (k = 0; k < nRepeatCount; k++)
690  {
691  CKNComplex temp = lpResult->pEigenVectorsForAMatrix[i].GetAt(k);
692  temp = temp + V.GetAt(k) * lpResult->pEigenVectors[i][nIterationIndex - 1];
693  lpResult->pEigenVectorsForAMatrix[i].SetAt(k, temp);
694  }
695  }
696 }
697 
708 bool CKNLanczosMethod::CheckAndDoSelectiveReorthogonalization(int nIterationCount, double *pAlpha, double *pBeta, double *pWj, double *pWjm1, double *pWjp1, double fANorm)
709 {
710  bool bDoSelectiveReorthogonalization = false;
711 
713  return bDoSelectiveReorthogonalization;
714 }
715 
719 void CKNLanczosMethod::DoSelectiveReorthogonalization(unsigned int nIterationCount)
720 {
722  return;
723 }
724 
725 void thomas_alg(double **T, double *initguess, double *app_evc, int iter)
726 { //[completed]
727 
728  int i;
729  double temp, tbeta = T[0][0];
730 
731  double *P = (double *)malloc(sizeof(double)*iter);
732 
733  //Forward Elimination
734  for (i = 0; i<iter - 1; i++)
735  {
736  if (fabs(T[1][i])>fabs(tbeta))
737  {
738  initguess[i + 1] -= initguess[i] * tbeta / T[1][i];
739  T[1][i + 1] -= T[0][i] * tbeta / T[1][i];
740  P[i] = 0;
741  tbeta = T[0][i + 1];
742  }
743  else
744  {
745  temp = initguess[i + 1];
746  initguess[i + 1] = initguess[i];
747  initguess[i] = temp;
748 
749  temp = T[1][i + 1];
750  T[1][i + 1] = T[0][i];
751  T[0][i] = temp;
752 
753  initguess[i + 1] -= initguess[i] * T[1][i] / tbeta;
754 
755  T[1][i + 1] -= T[0][i] * T[1][i] / tbeta;
756 
757  if (i != iter - 2)
758  {
759  P[i] = T[0][i + 1];
760  T[0][i + 1] = -P[i] * T[1][i] / tbeta;
761  }
762 
763  T[1][i] = tbeta;
764 
765  if (i != iter - 2)
766  tbeta = P[i];
767  }
768  }
769 
770  //Backward Substitution
771  for (i = iter - 1; i>-1; i--)
772  {
773  if (i == iter - 1)
774  app_evc[i] = initguess[i] / T[1][i];
775  else if (i == iter - 2)
776  {
777  app_evc[i] = initguess[i] - app_evc[i + 1] * T[0][i];
778  app_evc[i] = app_evc[i] / T[1][i];
779  }
780  else
781  {
782  app_evc[i] = initguess[i] - app_evc[i + 1] * T[0][i] - app_evc[i + 2] * P[i];
783  app_evc[i] = app_evc[i] / T[1][i];
784  }
785  }
786 
787  free(P);
788 }
789 
790 void inverse_iter(double *alpha, double *beta, double *app_evc, int iter, double app_eva)
791 { //[completed]
792  int i, j, k, iteration = 20;
793  double temp = 0, err, etol = 1e-13;
794 
795  double *initguess = (double *)malloc(sizeof(double)*iter); // Initial Guess
796  double **T = (double **)malloc(sizeof(double)* 2);
797 
798  for (i = 0; i<2; i++)
799  T[i] = (double *)malloc(sizeof(double)*iter);
800 
801  for (i = 0; i<iter; i++)
802  initguess[i] = 1 / sqrt((double)iter);
803 
804  for (k = 0; k<iteration; k++)
805  {
806  for (i = 0; i<2; i++)
807  {
808  if (i == 1)
809  {
810  for (j = 0; j<iter; j++)
811  T[i][j] = alpha[j] - app_eva;
812 
813  }
814  else
815  {
816  for (j = 0; j<iter; j++)
817  T[i][j] = beta[j];
818  }
819  }
820 
821  thomas_alg(T, initguess, app_evc, iter); // Solve (T-app_eva*I)*app_evc=initguess
822 
823  temp = 0;
824  for (i = 0; i<iter; i++)
825  temp += app_evc[i] * app_evc[i];
826 
827  for (i = 0; i<iter; i++)
828  app_evc[i] = app_evc[i] / sqrt(temp);
829 
830  temp = 0;
831  for (i = 0; i<iter; i++)
832  {
833  if (i == 0)
834  temp += app_evc[i] * (alpha[i] * app_evc[i] + beta[i] * app_evc[i + 1]);
835  else if (i == iter - 1)
836  temp += app_evc[i] * (beta[i - 1] * app_evc[i - 1] + alpha[i] * app_evc[i]);
837  else
838  temp += app_evc[i] * (beta[i - 1] * app_evc[i - 1] + alpha[i] * app_evc[i] + beta[i] * app_evc[i + 1]);
839  }
840 
841  err = app_eva - temp;
842  if (fabs(err) < etol)
843  break;
844 
845  if (k<iteration - 1)
846  {
847  for (i = 0; i<iter; i++)
848  initguess[i] = app_evc[i];
849  }
850  } // Iteration k
851 
852  for (i = 0; i<2; i++)
853  free(T[i]);
854  free(T);
855  free(initguess);
856 }
857 
858 #define SET_RESULT_TO_PARAMETER(pResultValue, pResultVector, nResultCount) \
859  pCalcuResult_Value = pResultValue; \
860  pCalcuResult_Vector = pResultVector; \
861  nCalculatedEigenValueCount = nResultCount
862 
872 bool CKNLanczosMethod::DoEigenValueSolving(int nIterationCount, double *pAlpha, double *pBeta, double fANorm, LPEIGENVALUE_RESULT lpResult, bool bFinal)
873 {
874  double *pCalcuResult_Value = NULL;
875  double *pCalcuResult_Vector = NULL;
876  int nConvergedEigenvalueCount = 0;
877  int nRangecheckedEigenvalueCount = 0;
878  int nNonSpuriousRitzValueCount = 0;
879  int nNonClustersValueCount = 0;
880  unsigned int nCalculatedEigenValueCount = 0, nCalculatedEigenValueCountBeforeConvergenceCheck;
881  bool bRtn = false;
882  bool *pValidEigenValue = NULL;
883  unsigned int i;
884 
885  if (NULL == lpResult)
886  throw ERROR_MALLOC;
887 
888  if (!CKNMPIManager::IsRootRank())
889  return false;
890 
891  nCalculatedEigenValueCount = EigenValueSolver(nIterationCount, pAlpha, pBeta, m_pEigenValues, m_pEigenVectors);
892  pCalcuResult_Value = m_pEigenValues;
893  pCalcuResult_Vector = m_pEigenVectors;
894  nCalculatedEigenValueCountBeforeConvergenceCheck = nCalculatedEigenValueCount;
895 
897  m_pEigenVectors = (double*)malloc(sizeof(double)*nCalculatedEigenValueCount*nIterationCount);
898  for (i = 0; i < nCalculatedEigenValueCount; ++i)
899  {
900  double *pEigenVector = (double*)malloc(sizeof(double)*nIterationCount);
901 
902  if (NULL == pEigenVector)
903  continue;
904 
905  inverse_iter(pAlpha + 1, pBeta + 2, pEigenVector, nIterationCount, m_pEigenValues[i]);
906  memcpy(m_pEigenVectors + nIterationCount * i, pEigenVector, sizeof(double)*nIterationCount);
907 
908  FREE_MEM(pEigenVector);
909  }
910 
911  pValidEigenValue = (bool*)malloc(sizeof(bool)*nCalculatedEigenValueCount);
912  for (i = 0; i < nCalculatedEigenValueCount; ++i)
913  pValidEigenValue[i] = true;
914 
915  if (DO_NOT_CONVERGENCE_CHECKING != m_fConvergenceTolerance)
916  {
917  nCalculatedEigenValueCount = ConvergenceCheckingEx(nCalculatedEigenValueCount, m_pEigenValues, m_pEigenVectors,
918  pValidEigenValue, fANorm, pBeta, nIterationCount);
919  }
920 
921  {
922  nCalculatedEigenValueCount = DistinguishClusterOfEigenvalueEx(nCalculatedEigenValueCount, m_pEigenValues, m_pEigenVectors,
923  pValidEigenValue, nIterationCount);
924  }
925 
926  if (nCalculatedEigenValueCount > 0)
927  {
928  IntegrateEigenvaluesEx(nIterationCount, lpResult, nCalculatedEigenValueCount, nCalculatedEigenValueCountBeforeConvergenceCheck, m_pEigenValues, m_pEigenVectors, pValidEigenValue);
929  }
930 
931  FREE_MEM(pValidEigenValue);
932 
933  if (lpResult->nEigenValueCount >= m_nEigenValueCount || true == bFinal)
934  bRtn = true;
935 
936  return bRtn;
937 }
938 
948 void CKNLanczosMethod::IntegrateEigenvaluesEx(int nIterationCount, LPEIGENVALUE_RESULT lpResult, unsigned int nCalculatedEigenValueCount, unsigned int nCalculatedEigenValueCountBeforeConvergenceCheck, double *pCalcuResult_Value, double *pCalcuResult_Vector, bool *pbValidEigenValue)
949 {
950  unsigned int i, j;
951  bool bSame = false;
952 
953  lpResult->nEigenVectorSize = nIterationCount;
954  if (0 == lpResult->nEigenValueCount)
955  {
956  //Initialize lpReulst members
957  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
958  lpResult->pEigenValues = (double*)malloc(sizeof(double)*nCalculatedEigenValueCount);
959  lpResult->pEigenVectors = (double**)malloc(sizeof(double*)*nCalculatedEigenValueCount);
960 
961  for (i = 0; i < nCalculatedEigenValueCount; i++)
962  {
963  lpResult->pEigenVectors[i] = (double*)malloc(sizeof(double)*nIterationCount);
964  }
965 
966  lpResult->pEigenValueFoundIteration = (unsigned int*)malloc(sizeof(unsigned int)*nCalculatedEigenValueCount);
967  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
968  }
969  else
970  {
971  int nSameCount = 0;
972  for (i = 0; i < nCalculatedEigenValueCountBeforeConvergenceCheck; i++)
973  {
974  if (false == pbValidEigenValue[i])
975  continue;
976  for (j = 0; j < lpResult->nEigenValueCount; j++)
977  {
978  if (CKNMatrixOperation::IsSame(pCalcuResult_Value[i], lpResult->pEigenValues[j], GENERAL_TOLERANCE))
979  {
980  nSameCount++;
981  continue;
982  }
983  }
984  }
985 
986  if (0 == nCalculatedEigenValueCount - nSameCount)
987  {
988  for (i = 0; i < lpResult->nEigenValueCount ; i++)
989  lpResult->pEigenVectors[i] = (double*)realloc(lpResult->pEigenVectors[i], sizeof(double)*nIterationCount);
990  }
991  else
992  {
993  int nAddSize = nCalculatedEigenValueCount - nSameCount;
994  lpResult->pEigenValues = (double*)realloc(lpResult->pEigenValues, sizeof(double)*(lpResult->nEigenValueCount + nAddSize));
995  lpResult->pEigenVectors = (double**)realloc(lpResult->pEigenVectors, sizeof(double*)*(lpResult->nEigenValueCount + nAddSize));
996  for (i = lpResult->nEigenValueCount; i < (lpResult->nEigenValueCount + nAddSize); ++i)
997  lpResult->pEigenVectors[i] = NULL;
998  for (i = 0; i < (lpResult->nEigenValueCount + nAddSize); i++)
999  lpResult->pEigenVectors[i] = (double*)realloc(lpResult->pEigenVectors[i], sizeof(double)*nIterationCount);
1000  lpResult->pEigenValueFoundIteration = (unsigned int*)realloc(lpResult->pEigenValueFoundIteration, sizeof(unsigned int)*(lpResult->nEigenValueCount + nAddSize));
1001  }
1002  }
1003 
1004  for (i = 0; i < nCalculatedEigenValueCountBeforeConvergenceCheck; i++)
1005  {
1006  if (false == pbValidEigenValue[i])
1007  continue;
1008 
1009  bSame = false;
1010  for (j = 0; j < lpResult->nEigenValueCount; j++)
1011  {
1012  if (CKNMatrixOperation::IsSame(pCalcuResult_Value[i], lpResult->pEigenValues[j], GENERAL_TOLERANCE))
1013  {
1014  bSame = true;
1015  break;
1016  }
1017  }
1018 
1019  if (!bSame)
1020  {
1021  lpResult->pEigenValues[lpResult->nEigenValueCount] = pCalcuResult_Value[i];
1022  lpResult->pEigenValueFoundIteration[lpResult->nEigenValueCount] = nIterationCount;
1023 #ifdef _WIN32
1024  lpResult->nMaxEigenValueFoundIteration = max((int)lpResult->nMaxEigenValueFoundIteration, (int)nIterationCount);
1025 #else //_WIN32
1026  lpResult->nMaxEigenValueFoundIteration = std::max((int)lpResult->nMaxEigenValueFoundIteration, (int)nIterationCount);
1027 #endif//
1028  memcpy(lpResult->pEigenVectors[lpResult->nEigenValueCount], pCalcuResult_Vector + (i * nIterationCount), sizeof(double)* nIterationCount);
1029  lpResult->nEigenValueCount++;
1030  }
1031  else
1032  {
1033  lpResult->pEigenValueFoundIteration[j] = nIterationCount;
1034  memcpy(lpResult->pEigenVectors[j], pCalcuResult_Vector + (i * nIterationCount), sizeof(double)* nIterationCount);
1035 #ifdef _WIN32
1036  lpResult->nMaxEigenValueFoundIteration = max((int)lpResult->nMaxEigenValueFoundIteration, (int)nIterationCount);
1037 #else //_WIN32
1038  lpResult->nMaxEigenValueFoundIteration = std::max((int)lpResult->nMaxEigenValueFoundIteration, (int)nIterationCount);
1039 #endif //_WIN32
1040  }
1041  }
1042 }
1043 
1051 void CKNLanczosMethod::IntegrateEigenvalues(int nIterationCount, LPEIGENVALUE_RESULT lpResult, unsigned int nCalculatedEigenValueCount, double *pCalcuResult_Value, double *pCalcuResult_Vector)
1052 {
1053  unsigned int i, j;
1054  bool bSame = false;
1055 
1056  for (i = 0; i < nCalculatedEigenValueCount; i++)
1057  {
1058  bSame = false;
1059  for (j = 0; j < lpResult->nEigenValueCount; j++)
1060  {
1061  if (CKNMatrixOperation::IsSame(pCalcuResult_Value[i], lpResult->pEigenValues[j], GENERAL_TOLERANCE))
1062  {
1063  bSame = true;
1064  break;
1065  }
1066  }
1067 
1068  if (!bSame)
1069  {
1070  if (lpResult->nEigenValueCountForMemeory == lpResult->nEigenValueCount)
1071  {
1072  lpResult->pEigenValues = (double*)realloc(lpResult->pEigenValues, sizeof(double)*(lpResult->nEigenValueCount * 2));
1073  lpResult->nEigenValueCountForMemeory = lpResult->nEigenValueCount * 2;
1074  lpResult->pEigenValueFoundIteration = (unsigned int*)realloc(lpResult->pEigenValueFoundIteration, sizeof(unsigned int)*(lpResult->nEigenValueCount * 2));
1075  //lpResult->pEigenVectors = (double*)realloc(lpResult->pEigenVectors, sizeof(double)*lpResult->nEigenValueCount * 2 * m_nIterationCount);
1076  lpResult->nEigenValueCountForMemeory = lpResult->nEigenValueCount * 2;
1077  }
1078  lpResult->pEigenValues[lpResult->nEigenValueCount] = pCalcuResult_Value[i];
1079  lpResult->pEigenValueFoundIteration[lpResult->nEigenValueCount] = nIterationCount;
1080 #ifdef _WIN32
1081  lpResult->nMaxEigenValueFoundIteration = max((int)lpResult->nMaxEigenValueFoundIteration, (int)nIterationCount);
1082 #else //_WIN32
1083  lpResult->nMaxEigenValueFoundIteration = std::max((int)lpResult->nMaxEigenValueFoundIteration, (int)nIterationCount);
1084 #endif//
1085  memcpy(lpResult->pEigenVectors + (lpResult->nEigenValueCount*m_nIterationCount), pCalcuResult_Vector + (i * nIterationCount), sizeof(double)* nIterationCount);
1086  lpResult->nEigenValueCount++;
1087  }
1088  }
1089 }
1090 
1101 int CKNLanczosMethod::ConvergenceCheckingEx(int nEigenValueCount, double *pEigenValues, double *pEiegnVectors, bool *pbValidEigenValue, double fANorm, double *pBeta, int nIterationCount)
1102 {
1103  if (0 == nEigenValueCount)
1104  return 0;
1105 
1106  double fTotal = m_fConvergenceTolerance;
1107  double fResidual;
1108  int i, nConvergedCount = nEigenValueCount;
1109 
1110  for (i = 0; i < nEigenValueCount; i++)
1111  {
1112  if (false == pbValidEigenValue[i])
1113  continue;
1114 
1115  fResidual = fabs(pBeta[nIterationCount + 1] * pEiegnVectors[i*nIterationCount + nIterationCount - 1]);
1116 
1117  if (fResidual >= fTotal)
1118  {
1119  pbValidEigenValue[i] = false;
1120  nConvergedCount--;
1121  }
1122  }
1123 
1124  return nConvergedCount;
1125 }
1126 
1138 int CKNLanczosMethod::ConvergenceChecking(int nEigenValueCount, double *pEigenValues, double *pEiegnVectors, double *pConvergedEigenValues, double *pConvergedEigenVectors, double fANorm, double *pBeta, int nIterationCount)
1139 {
1140  if (0 == nEigenValueCount)
1141  return 0;
1142 
1143  if (NULL == pConvergedEigenValues || NULL == pConvergedEigenVectors || NULL == m_pConvergedIndex)
1144  throw ERROR_MALLOC;
1145 
1146  double fTotal = m_fConvergenceTolerance;
1147  double fResidual;
1148  int i, nConvergedCount = 0;
1149 
1150  for (i = 0; i < nEigenValueCount; i++)
1151  {
1152  fResidual = fabs(pBeta[nIterationCount + 1] * pEiegnVectors[i*nIterationCount + nIterationCount - 1]);
1153 
1154  if (fResidual < fTotal)
1155  m_pConvergedIndex[nConvergedCount++] = i;
1156  }
1157 
1158  ExtractDoubleValues(pConvergedEigenValues, pEigenValues, nEigenValueCount, m_pConvergedIndex, nConvergedCount, false);
1159  ExtractDoubleVector(nIterationCount, pConvergedEigenVectors, pEiegnVectors, nEigenValueCount, m_pConvergedIndex, nConvergedCount, false);
1160 
1161  return nConvergedCount;
1162 }
1163 
1173 int CKNLanczosMethod::RangeChecking(int nEigenValueCount, double *pEigenValues, double *pEiegnVectors, double *pRangeCheckingEigenValues, double *pRangeCheckingVectors, int nIterationCount)
1174 {
1175  int i, nRangecheckedCount = 0;
1176 
1177  if (NULL == pRangeCheckingEigenValues || NULL == pRangeCheckingVectors || NULL == m_pRangecheckedIndex)
1178  throw ERROR_MALLOC;
1179 
1180  for (i = 0; i < nEigenValueCount; i++)
1181  {
1182  if (pEigenValues[i] >= m_fEigenvalueMin && pEigenValues[i] <= m_fEignevalueMax)
1183  m_pRangecheckedIndex[nRangecheckedCount++] = i;
1184  }
1185 
1186  ExtractDoubleValues(pRangeCheckingEigenValues, pEigenValues, nEigenValueCount, m_pRangecheckedIndex, nRangecheckedCount, false);
1187  ExtractDoubleVector(nIterationCount, pRangeCheckingVectors, pEiegnVectors, nEigenValueCount, m_pRangecheckedIndex, nRangecheckedCount, false);
1188 
1189  return nRangecheckedCount;
1190 }
1191 
1202 int CKNLanczosMethod::SpuriousRitzValueChecking(int nEigenValueCount, double *pEigenValues, double *pEigenVectors, double *pNonSpuriousValues, double *pNonSpuriousVectors, double fANorm, int nIterationCount)
1203 {
1204  if (0 == nEigenValueCount)
1205  return 0;
1206 
1207  if (NULL == pNonSpuriousValues || NULL == pNonSpuriousVectors || NULL == m_pNonSpuriousValueIndex)
1208  throw ERROR_MALLOC;
1209 
1210  double eps = 1e-8;
1211  double fTotal = eps;
1212  int i, nNonSpuriousValue = 0;
1213 
1214  for (i = 0; i < nEigenValueCount; i++)
1215  {
1216  if (fabs(pEigenVectors[i*nIterationCount]) > fTotal)
1217  m_pNonSpuriousValueIndex[nNonSpuriousValue++] = i;
1218  }
1219 
1220  ExtractDoubleValues(pNonSpuriousValues, pEigenValues, nEigenValueCount, m_pNonSpuriousValueIndex, nNonSpuriousValue, false);
1221  ExtractDoubleVector(nIterationCount, pNonSpuriousVectors, pEigenVectors, nEigenValueCount, m_pNonSpuriousValueIndex, nNonSpuriousValue, false);
1222 
1223  return nNonSpuriousValue;
1224 }
1233 int CKNLanczosMethod::DistinguishClusterOfEigenvalueEx(int nEigenValueCount, double *pEigenValues, double *pEigenVectors, bool *pbValidEigenValues, int nIterationCount)
1234 {
1235  double eps = 1e-8;
1236  int i, j, nNonClustersValue = nEigenValueCount;
1237 
1238  if (0 == nEigenValueCount)
1239  return 0;
1240 
1241  for (i = 0; i < nEigenValueCount - 1; i++)
1242  {
1243  if (!pbValidEigenValues[i])
1244  continue;
1245 
1246  for (j = i + 1; j < nEigenValueCount; j++)
1247  {
1248  if (!pbValidEigenValues[j])
1249  continue;
1250 
1251  if (fabs(pEigenValues[i] - pEigenValues[j]) < eps)
1252  {
1253  pbValidEigenValues[j] = false;
1254  nNonClustersValue--;
1255  }
1256  }
1257  }
1258 
1259  return nNonClustersValue;
1260 }
1261 
1271 int CKNLanczosMethod::DistinguishClusterOfEigenvalue(int nEigenValueCount, double *pEigenValues, double *pEigenVectors, double *pNonClustersValues, double *pNonClustersVectors, int nIterationCount)
1272 {
1273  double eps = 1e-8;
1274  int i, j, nNonClustersValue = 0;
1275 
1276  if (0 == nEigenValueCount)
1277  return 0;
1278 
1279  if (NULL == pNonClustersValues || NULL == m_pCheckNonClusterValue || NULL == m_pNonClustersValueIndex)
1280  throw ERROR_MALLOC;
1281 
1282  for (i = 0; i < nEigenValueCount; i++)
1283  m_pCheckNonClusterValue[i] = true;
1284 
1285  for (i = 0; i < nEigenValueCount - 1; i++)
1286  {
1287  if (!m_pCheckNonClusterValue[i])
1288  continue;
1289 
1290  for (j = i + 1; j < nEigenValueCount; j++)
1291  {
1292  if (fabs(pEigenValues[i] - pEigenValues[j]) < eps)
1293  {
1294  m_pCheckNonClusterValue[j] = false;
1295  }
1296  }
1297  }
1298 
1299  for (i = 0; i < nEigenValueCount; i++)
1300  {
1301  if (m_pCheckNonClusterValue[i])
1302  m_pNonClustersValueIndex[nNonClustersValue++] = i;
1303  }
1304 
1305  ExtractDoubleValues(pNonClustersValues, pEigenValues, nEigenValueCount, m_pNonClustersValueIndex, nNonClustersValue, false);
1306  if (NULL != pNonClustersVectors)
1307  ExtractDoubleVector(nIterationCount, pNonClustersVectors, pEigenVectors, nEigenValueCount, m_pNonClustersValueIndex, nNonClustersValue, false);
1308 
1309  return nNonClustersValue;
1310 }
1311 
1312 
1321 int CKNLanczosMethod::EigenValueSolver(unsigned int nIterationCount, double *pAlpha, double *pBeta, double *pEigenValues, double *pEigenVectors)
1322 {
1323  MKL_INT n, lda, ldz, il, iu, lwork;
1324  MKL_INT nFoundEigenValueCount;
1325  MKL_INT info;
1326  double vl, vu;
1327  double abstol = 1.0e-8;
1328  double *work;
1329  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::EVALUE_MALLOC);
1330  //MKL_INT *iwork = (MKL_INT*)malloc(sizeof(MKL_INT)* 5 * nIterationCount);
1331  MKL_INT *iwork = (MKL_INT*)malloc(sizeof(MKL_INT)* 10 * nIterationCount);
1332  MKL_INT *ifail = (MKL_INT*)malloc(sizeof(MKL_INT)* nIterationCount);
1334  int nsplit = 0; /* number of diagonal blocks in matrix */
1335  int *iblock = (int *)malloc(sizeof(int)*nIterationCount);
1336  int *isplit = (int *)malloc(sizeof(int)*nIterationCount);
1337 
1338  //work = (double*)malloc(sizeof(double)*nIterationCount*4);
1339  work = (double*)malloc(sizeof(double)*nIterationCount * 10);
1340  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::EVALUE_MALLOC);
1341  n = ldz = lda = nIterationCount;
1342  il = 1;
1343  iu = nIterationCount;
1344  vl = m_fEigenvalueMin;
1345  vu = m_fEignevalueMax;
1346 
1348  lwork = -1;
1349  dstebz("V", "E", &n, &vl, &vu, &il, &iu, &abstol, pAlpha + 1, pBeta + 2, &nFoundEigenValueCount, &nsplit, pEigenValues, iblock, isplit, work, iwork, &info);
1350 
1351  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::EVALUE_FREE_MEM);
1352  FREE_MEM(iwork);
1353  FREE_MEM(work);
1354  FREE_MEM(ifail);
1355  FREE_MEM(iblock);
1356  FREE_MEM(isplit);
1357  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::EVALUE_FREE_MEM);
1358 
1359  return nFoundEigenValueCount;
1360 }
1361 
1367 double* CKNLanczosMethod::BuildTMatrix(unsigned int nOrder, double *pAlpha, double *pBeta)
1368 {
1369  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
1370  double *pTMatrix = (double*)malloc(sizeof(double)*nOrder*nOrder);
1371  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
1372 
1373  unsigned int k;
1374 
1375  if (NULL == pTMatrix)
1376  {
1377  throw ERROR_MALLOC;
1378  return NULL;
1379  }
1380 
1381  memset(pTMatrix, 0, sizeof(double)*nOrder*nOrder);
1382 
1383  for (k = 1; k <= nOrder; k++)
1384  {
1385  pTMatrix[((k - 1) * nOrder) + (k - 1)] = pAlpha[k];
1386  if (k != nOrder)
1387  pTMatrix[(k - 1)*nOrder + k] = pBeta[k + 1];
1388  if (k != 1)
1389  pTMatrix[(k - 1)*nOrder + (k - 2)] = pBeta[k];
1390  }
1391 
1392  return pTMatrix;
1393 }
1394 
1395 void CKNLanczosMethod::InitVariables()
1396 {
1397  if (NULL != m_pV)
1398  {
1399  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::FREE_MEM);
1400  delete m_pV;
1401  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::FREE_MEM);
1402  m_pV = NULL;
1403  }
1404 
1405  m_pAMatrix = NULL;
1406  m_nMatrixSize = 0;
1407  m_nIterationCount = 0;
1408  m_nEigenValueCheckInterval = 0;
1409  m_nEigenValueCount = 0;
1410  m_bReorthogonalization = false;
1411 }
1412 
1421 void CKNLanczosMethod::ExtractDoubleValues(double *pTarget, double *pSource, unsigned int nSrcCount, int *pFilter, unsigned int nFilterCount, bool bExclusive)
1422 {
1423  unsigned int i;
1424 
1425  if (bExclusive)
1426  {
1427  int nSrcIndex = 0;
1428  int nTargetIndex = 0;
1429  for (i = 0; i < nSrcCount; i++)
1430  {
1431  if (i != pFilter[nSrcIndex])
1432  pTarget[nTargetIndex++] = pSource[i];
1433  else
1434  nSrcIndex++;
1435  }
1436  }
1437  else
1438  {
1439  for (i = 0; i < nFilterCount; i++)
1440  pTarget[i] = pSource[pFilter[i]];
1441  }
1442 }
1443 
1453 void CKNLanczosMethod::ExtractDoubleVector(unsigned int nVectorsize, double *pTarget, double *pSource, unsigned int nSrcCount, int *pFilter, unsigned int nFilterCount, bool bExclusive)
1454 {
1455  unsigned int i;
1456 
1457  if (bExclusive)
1458  {
1459  int nSrcIndex = 0;
1460  int nTargetIndex = 0;
1461  for (i = 0; i < nSrcCount; i++)
1462  {
1463  if (i != pFilter[nSrcIndex])
1464  memcpy(pTarget + (nTargetIndex*nVectorsize), pSource + (i * nVectorsize), sizeof(double)* nVectorsize);
1465  else
1466  nSrcIndex++;
1467  }
1468  }
1469  else
1470  {
1471  for (i = 0; i < nFilterCount; i++)
1472  memcpy(pTarget + (i*nVectorsize), pSource + (pFilter[i] * nVectorsize), sizeof(double)* nVectorsize);
1473  }
1474 }
1475 
1480 void CKNLanczosMethod::ReleaseResult(LPEIGENVALUE_RESULT lpResult, bool bReleaseStruct)
1481 {
1482  unsigned int i;
1483  if (NULL == lpResult)
1484  return;
1485 
1486  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::FREE_MEM);
1487 
1488  FREE_MEM(lpResult->pEigenValues);
1489  FREE_MEM(lpResult->pEigenValueFoundIteration);
1490 
1491  if (NULL != lpResult->pWaveFunctions)
1492  {
1493  for (i = 0; i < lpResult->nEigenValueCount / 10; i++)
1494  lpResult->pWaveFunctions[i].Finalize();
1495 
1496  delete[] lpResult->pWaveFunctions;
1497  lpResult->pWaveFunctions = NULL;
1498  }
1499 
1500  if (NULL != lpResult->pEigenVectorsForAMatrix)
1501  {
1502  for (i = 0; i < lpResult->nEigenValueCount; i++)
1503  lpResult->pEigenVectorsForAMatrix[i].Finalize();
1504 
1505  delete[] lpResult->pEigenVectorsForAMatrix;
1506  lpResult->pEigenVectorsForAMatrix = NULL;
1507  }
1508 
1509  if( NULL != lpResult->pEigenVectors)
1510  {
1511  for (i = 0; i < lpResult->nEigenValueCount - lpResult->nDegeneratedEigenValueCount ; i++)
1512  FREE_MEM(lpResult->pEigenVectors[i]);
1513 
1514  FREE_MEM(lpResult->pEigenVectors);
1515  }
1516 
1517  if (bReleaseStruct)
1518  FREE_MEM(lpResult);
1519  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::FREE_MEM);
1520 }
1521 
1522 void CKNLanczosMethod::StopIteration()
1523 {
1524  CKNLanczosMethod::m_bStop = true;
1525 }
1526 
1530 void CKNLanczosMethod::BuildWaveFunction(LPEIGENVALUE_RESULT lpResult)
1531 {
1532  if (0 != m_nMatrixSize % 10 || 0 != CKNMPIManager::GetCurrentLoadBalanceCount() % 10)
1533  return;
1534 
1535  unsigned int i, j, nIndex = 0;
1536  CKNComplex tempResult, complexNumber;
1537  double fTempResult;
1538 
1539  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MALLOC);
1540  lpResult->pWaveFunctions = new CKNMatrixOperation::CKNVector[lpResult->nEigenValueCount];
1541  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MALLOC);
1542  for (i = 0; i < lpResult->nEigenValueCount; i++)
1543  {
1544  nIndex = 0;
1545  fTempResult = 0.;
1546  lpResult->pWaveFunctions[i].SetSize(CKNMPIManager::GetCurrentLoadBalanceCount() / 10);
1547  for (j = 0; j < CKNMPIManager::GetCurrentLoadBalanceCount(); j++)
1548  {
1549  complexNumber = lpResult->pEigenVectorsForAMatrix[i].GetAt(j);
1550  double absoluteValue = complexNumber.GetAbsolute();
1551  fTempResult += (absoluteValue*absoluteValue);
1552  if (9 == j % 10)
1553  {
1554  lpResult->pWaveFunctions[i].SetAt(nIndex++, fTempResult, 0);
1555  fTempResult = 0.;
1556  }
1557  }
1558  }
1559 }
1560 
1565 void CKNLanczosMethod::DoResidualCheck(CKNMatrixOperation::CKNCSR *pAMatrix, LPEIGENVALUE_RESULT lpResult)
1566 {
1567  int i, j;
1568  CKNMatrixOperation::CKNVector vectorResult1, vectorResult2;
1569  CKNComplex result;
1570  double fEigenValue;
1571  double fNorm = 0.0;
1572  std::vector<bool> vectorResidualCheck;
1573  bool bFoundNoAnswer = false;
1574  std::vector<int> vectorResidualAnswer;
1575  int nResidualAnswerCount = 0;
1576 
1577  SHOW_SIMPLE_MSG("-Residual Checking...\n");
1578 
1579 #ifdef DISABLE_MPI_ROUTINE
1580  vectorResult1.SetSize(pAMatrix->GetColumnCount());
1581  vectorResult2.SetSize(pAMatrix->GetColumnCount());
1582 #else //DISABLE_MPI_ROUTINE
1583  vectorResult1.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
1584  vectorResult2.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
1585 #endif //DISABLE_MPI_ROUTINE
1586  for( i = 0 ; i < lpResult->nEigenValueCount ; ++i )
1587  {
1588  CKNMatrixOperation::MVMul(pAMatrix, &lpResult->pEigenVectorsForAMatrix[i], &vectorResult1);
1589  vectorResult2 = lpResult->pEigenVectorsForAMatrix[i];
1590  if( CKNMPIManager::IsRootRank() )
1591  fEigenValue = lpResult->pEigenValues[i];
1592  CKNMPIManager::BroadcastDouble(&fEigenValue, 1);
1593  vectorResult2.ScalarMultiple(fEigenValue);
1594  //CKNMatrixOperation::VVDot(&vectorResult1, &vectorResult2, &result);
1595  vectorResult1.MinusVector(&vectorResult2);
1596 #ifdef DISABLE_MPI_ROUTINE
1597  fNorm = vectorResult1.GetNorm();
1598 #else //DISABLE_MPI_ROUTINE
1599  fNorm = vectorResult1.GetNorm(true);
1600 #endif DISABLE_MPI_ROUTINE
1601 
1602  if (CKNMatrixOperation::IsSame(fNorm, 0.0, GENERAL_TOLERANCE))
1603  vectorResidualCheck.push_back(true);
1604  else
1605  {
1606  vectorResidualCheck.push_back(false);
1607  bFoundNoAnswer = true;
1608  }
1609  }
1610 
1611  if (!bFoundNoAnswer)
1612  {
1613  vectorResult1.Finalize();
1614  vectorResult2.Finalize();
1615  return;
1616  }
1617 
1618  for( i = 0 ; i < lpResult->nEigenValueCount ; ++ i )
1619  {
1620  if( vectorResidualCheck[i] )
1621  {
1622  vectorResidualAnswer.push_back(i);
1623  nResidualAnswerCount++;
1624  }
1625  }
1626 
1627  for( i = 0 ; i < nResidualAnswerCount; ++ i )
1628  {
1629  if( i == vectorResidualAnswer[i] )
1630  continue;
1631 
1632  lpResult->pEigenValues[i] = lpResult->pEigenValues[vectorResidualAnswer[i]];
1633  lpResult->pEigenValueFoundIteration[i] = lpResult->pEigenValueFoundIteration[vectorResidualAnswer[i]];
1634  lpResult->pEigenVectorsForAMatrix[i] = lpResult->pEigenVectorsForAMatrix[vectorResidualAnswer[i]];
1635 
1636  if( NULL != lpResult->pWaveFunctions )
1637  lpResult->pWaveFunctions[i] = lpResult->pWaveFunctions[vectorResidualAnswer[i]];
1638  }
1639 
1640  lpResult->nEigenValueCount = nResidualAnswerCount;
1641 
1642  vectorResult1.Finalize();
1643  vectorResult2.Finalize();
1644 }
1645 
1654 void CKNLanczosMethod::SaveLanczosResult(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, bool bCalcuEigenvalue, bool bWaveFunction, double *pKValue, int nRepeatCount)
1655 {
1656  if (NULL == lpResult || NULL == pKValue)
1657  return;
1658  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::FILEIO);
1659  FILE *out;
1660  char szFileName[1024], szBuffer[1024];
1661  std::string writeString;
1662  unsigned int i, j, k;
1663  char szFileOpt[3] = "wt";
1664 
1665  if( 0 != nRepeatCount )
1666  strcpy(szFileOpt, "at");
1667 
1668 #ifdef _WIN32
1669  _mkdir("result");
1670 #else
1671  mkdir("result", 0777);
1672 #endif
1673 
1674  if( CKNMPIManager::IsDeflationRoot() && CKNMPIManager::IsRootRank())
1675  {
1676 #ifdef _WIN32
1677  if (NULL != (out = fopen("result\\eigenvalues.txt", szFileOpt)))
1678 #else _WIN32
1679  if (NULL != (out = fopen("result/eigenvalues.txt", szFileOpt)))
1680 #endif //_WIN32
1681  {
1682  int nEigenValueCount;
1683  if( NULL != pKValue )
1684  sprintf(szBuffer, "---[ev #%d %8.6f %8.6f %8.6f]---\n\n", nRepeatCount, pKValue[0], pKValue[1], pKValue[2]);
1685  else
1686  sprintf(szBuffer, "---[ev]---\n\n");
1687 
1688  fputs(szBuffer, out);
1689  for (i = 0; i < lpResult->nEigenValueCount; i++)
1690  {
1691 
1692  if( DEGENERATED_INDEX != lpResult->pEigenValueFoundIteration[i] )
1693  sprintf(szBuffer, "[ev %2d] %18.16f - Iteration [#%7d]\n", i, lpResult->pEigenValues[i], lpResult->pEigenValueFoundIteration[i]);
1694  else
1695  sprintf(szBuffer, "[ev %2d] %18.16f - Degenerated eigenvalue\n", i, lpResult->pEigenValues[i], lpResult->pEigenValueFoundIteration[i]);
1696 
1697  fputs(szBuffer, out);
1698  }
1699 
1700  /*nEigenValueCount = i;
1701  if( lpResult->nDegeneratedEigenValueCount > 0 )
1702  {
1703  for(i = 0; i < lpResult->nDegeneratedEigenValueCount; ++i)
1704  {
1705  sprintf(szBuffer, "[ev %2d] %18.16f - Degenerated eigenvalue\n", nEigenValueCount + i, lpResult->pDegeneratedEigenValues[i]);
1706  fputs(szBuffer, out);
1707  }
1708  }*/
1709 
1710  fputs("\n\n", out);
1711  fclose(out);
1712  }
1713  }
1714 
1715  if (bCalcuEigenvalue && NULL != lpResult->pEigenVectorsForAMatrix)
1716  {
1717  if( CKNMPIManager::IsDeflationRoot() )
1718  {
1719  for (j = 0; j < lpResult->nEigenValueCount; j++)
1720  {
1721 #ifdef _WIN32
1722  sprintf(szFileName, "result\\eigenvector_%02d_%02d.txt", nRepeatCount, j);
1723 #else _WIN32
1724  sprintf(szFileName, "result/eigenvector_%02d_%02d.txt", nRepeatCount, j);
1725 #endif //_WIN32
1726  for( k = 0; k < CKNMPIManager::GetTotalNodeCount() ; ++k)
1727  {
1728  if( k == CKNMPIManager::GetCurrentRank() )
1729  {
1730  if (NULL != (out = fopen(szFileName, "at")))
1731  {
1732 
1733  for (i = 0; i < lpResult->pEigenVectorsForAMatrix[j].GetSize(); i++)
1734  {
1735  sprintf(szBuffer, "%16.16f %16.16f\n",
1736  lpResult->pEigenVectorsForAMatrix[j].GetAt(i).GetRealNumber(),
1737  lpResult->pEigenVectorsForAMatrix[j].GetAt(i).GetImaginaryNumber());
1738 
1739  writeString += szBuffer;
1740 
1741  if (i % 100)
1742  {
1743  fputs(writeString.c_str(), out);
1744  writeString.clear();
1745  }
1746  }
1747 
1748  if (!writeString.empty())
1749  {
1750  fputs(writeString.c_str(), out);
1751  writeString.clear();
1752  }
1753 
1754  fclose(out);
1755  }
1756  }
1757 #ifndef DISABLE_MPI_ROUTINE
1758  CKNMPIManager::Barrier();
1759 #endif //DISABLE_MPI_ROUTINE
1760  }
1761  }
1762  }
1763  }
1764 
1765  if (bWaveFunction && NULL != lpResult->pWaveFunctions)
1766  {
1767  if( CKNMPIManager::IsDeflationRoot() )
1768  {
1769  for (j = 0; j < lpResult->nEigenValueCount; j++)
1770  {
1771 #ifdef _WIN32
1772  sprintf(szFileName, "result\\wavefunction_%02d_%02d.txt", nRepeatCount, j);
1773 #else _WIN32
1774  sprintf(szFileName, "result/wavefunction_%02d_%02d.txt", nRepeatCount, j);
1775 #endif //_WIN32
1776  for( k = 0; k < CKNMPIManager::GetTotalNodeCount() ; ++k)
1777  {
1778  if( k == CKNMPIManager::GetCurrentRank() )
1779  {
1780  if (NULL != (out = fopen(szFileName, "at")))
1781  {
1782  for (i = 0; i < lpResult->pWaveFunctions[j].GetSize(); i++)
1783  {
1784  sprintf(szBuffer, "%16.16f\n",
1785  lpResult->pWaveFunctions[j].GetAt(i).GetRealNumber());
1786 
1787  writeString += szBuffer;
1788 
1789  if (i % 100)
1790  {
1791  fputs(writeString.c_str(), out);
1792  writeString.clear();
1793  }
1794  }
1795 
1796  if (!writeString.empty())
1797  {
1798  fputs(writeString.c_str(), out);
1799  writeString.clear();
1800  }
1801  fclose(out);
1802  }
1803  }
1804 #ifndef DISABLE_MPI_ROUTINE
1805  CKNMPIManager::Barrier();
1806 #endif //DISABLE_MPI_ROUTINE
1807  }
1808  }
1809  }
1810  }
1811 
1812  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::FILEIO);
1813 }
1814 
1822 void CKNLanczosMethod::ShowLanczosResult(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, bool bCalculateEigenVectors, bool bCalculateWaveFunction, double *pKValue, int nRepeatCount)
1823 {
1824  char szMsg[1024];
1825  int i;
1826 
1827  if( !CKNMPIManager::IsDeflationRoot() || !CKNMPIManager::IsRootRank())
1828  return;
1829 
1830  double fEvalTime = CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::EVALUE) - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::EVALUE_MALLOC) - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::EVALUE_FREE_MEM);
1831  double fCommTime = CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::COMM) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MV_COMM) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::VV_COMM);
1832  double fMemTime = CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MALLOC) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::FREE_MEM) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::EVALUE_MALLOC) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::EVALUE_FREE_MEM) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MV_MALLOC) + CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MV_FREE_MEM) ;
1833  double fMVMulTime = CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MVMUL) - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MV_MALLOC) - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MV_FREE_MEM) - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::MV_COMM);
1834  double fVVTime = CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::VVDOT) - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::VV_COMM);
1835  double fComputingTime = CKNTimeMeasurement::GetTotalTakeTime() - fEvalTime - fCommTime - fMemTime - fMVMulTime - fVVTime - CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::FILEIO);
1836 
1837 #ifdef DOING_MEASUREMENT
1838  CKNIPCCUtility::ShowMsg("\n--------------------------------------\nTime evaluation\n\n");
1839  sprintf(szMsg, "%d nodes used\nTotal time [ %f\tsec ]\nComputing [ %lf\tsec ]\nEvalue takes [ %lf\tsec ]\nMPI takes [ %lf\tsec ]\nMVMul takes [ %lf\tsec ]\nVVDot takes [ %lf\tsec ]\nMem Op takes [ %lf\tsec ]\nResult written [ %lf\tsec ]\n",
1840  CKNMPIManager::GetTotalNodeCount(),
1841  CKNTimeMeasurement::GetTotalTakeTime(),
1842  fComputingTime,
1843  fEvalTime,
1844  fCommTime,
1845  fMVMulTime,
1846  fVVTime,
1847  fMemTime,
1848  CKNTimeMeasurement::GetTakeTime(CKNTimeMeasurement::FILEIO));
1849 
1850  CKNIPCCUtility::ShowMsg(szMsg);
1851  CKNIPCCUtility::ShowMsg("--------------------------------------\n");
1852 #endif //DOING_MEASUREMENT
1853 
1854  if (NULL != lpResult)
1855  {
1856  int nEigenValueCount;
1857  CKNIPCCUtility::ShowMsg("\n--------------------------------------------------\nEigen values\n\n");
1858  for (i = 0; i < (int)lpResult->nEigenValueCount; ++i)
1859  {
1860  if( DEGENERATED_INDEX != lpResult->pEigenValueFoundIteration[i] )
1861  sprintf(szMsg, "[ev %2d] %18.16f - Iteration [#%7d]\n", i, lpResult->pEigenValues[i], lpResult->pEigenValueFoundIteration[i]);
1862  else
1863  sprintf(szMsg, "[ev %2d] %18.16f - Degenerated eigenvalue\n", i, lpResult->pEigenValues[i]);
1864  CKNIPCCUtility::ShowMsg(szMsg);
1865  }
1866  CKNIPCCUtility::ShowMsg("--------------------------------------------------\n");
1867  }
1868 }
1869 
1873 void CKNLanczosMethod::RecalcuWaveFunction(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult)
1874 {
1875 
1876  if( !CKNMPIManager::IsDeflationRoot() )
1877  return;
1878 
1879  unsigned int i, j, nIndex = 0;
1880  CKNComplex tempResult, complexNumber;
1881  double fTempResult;
1882 
1883  for (i = 0; i < lpResult->nEigenValueCount; i++)
1884  {
1885  nIndex = 0;
1886  fTempResult = 0.;
1887  for (j = 0; j < CKNMPIManager::GetCurrentLoadBalanceCount(); j++)
1888  {
1889  complexNumber = lpResult->pEigenVectorsForAMatrix[i].GetAt(j);
1890  double absoluteValue = complexNumber.GetAbsolute();
1891  fTempResult += (absoluteValue*absoluteValue);
1892  if (9 == j % 10)
1893  {
1894  lpResult->pWaveFunctions[i].SetAt(nIndex++, fTempResult, 0);
1895  fTempResult = 0.;
1896  }
1897  }
1898  }
1899 
1900 }
1901 
1908 void CKNLanczosMethod::AppendEigenValue(LPEIGENVALUE_RESULT lpResult, double fEigenValue, unsigned int nFindIteration, bool bInsertFirst)
1909 {
1910  int i;
1911 
1912  lpResult->nEigenValueCount++;
1913  lpResult->pEigenValues = (double*)realloc(lpResult->pEigenValues, sizeof(double) * lpResult->nEigenValueCount);
1914  if (bInsertFirst)
1915  {
1916  for( i = lpResult->nEigenValueCount - 1 ; i > 0 ; --i )
1917  lpResult->pEigenValues[i] = lpResult->pEigenValues[i-1];
1918 
1919  lpResult->pEigenValues[0] = fEigenValue;
1920  }
1921  else
1922  lpResult->pEigenValues[lpResult->nEigenValueCount-1] = fEigenValue;
1923 
1924  lpResult->pEigenValueFoundIteration = (unsigned int*)realloc(lpResult->pEigenValueFoundIteration, sizeof(unsigned int)*lpResult->nEigenValueCount);
1925  if (bInsertFirst)
1926  {
1927  for( i = lpResult->nEigenValueCount - 1 ; i > 0 ; --i )
1928  lpResult->pEigenValueFoundIteration[i] = lpResult->pEigenValueFoundIteration[i-1];
1929 
1930  lpResult->pEigenValueFoundIteration[0] = nFindIteration;
1931  }
1932  else
1933  lpResult->pEigenValueFoundIteration[lpResult->nEigenValueCount-1] = nFindIteration;
1934 }
1935 
1941 void CKNLanczosMethod::AppendEigenVector(LPEIGENVALUE_RESULT lpResult, CKNMatrixOperation::CKNVector *pEigenVector, bool bInsertFirst)
1942 {
1943  unsigned int i, nNewIndex;
1944  CKNMatrixOperation::CKNVector *pEigenVectorsForAMatrix = new CKNMatrixOperation::CKNVector[lpResult->nEigenValueCount];
1945  CKNMatrixOperation::CKNVector *pWaveFunctions = NULL;
1946  unsigned int nAdjust = 0;
1947 
1948 
1949  if (bInsertFirst)
1950  nAdjust = 1;
1951 
1952 
1953  if( NULL != lpResult->pWaveFunctions )
1954  pWaveFunctions = new CKNMatrixOperation::CKNVector[lpResult->nEigenValueCount];
1955 
1956  for( i = 0 + nAdjust; i < lpResult->nEigenValueCount - 1 + nAdjust; ++ i)
1957  {
1958  pEigenVectorsForAMatrix[i].SetSize(lpResult->pEigenVectorsForAMatrix[i-nAdjust].GetSize());
1959  pEigenVectorsForAMatrix[i] = lpResult->pEigenVectorsForAMatrix[i-nAdjust];
1960 
1961  if( NULL != lpResult->pWaveFunctions )
1962  {
1963  pWaveFunctions[i].SetSize(lpResult->pWaveFunctions[i-nAdjust].GetSize());
1964  pWaveFunctions[i] = lpResult->pWaveFunctions[i-nAdjust];
1965 
1966  if (bInsertFirst)
1967  pWaveFunctions[0].SetSize(pEigenVector->GetSize()/10);
1968  else
1969  pWaveFunctions[lpResult->nEigenValueCount-1].SetSize(pEigenVector->GetSize()/10);
1970  }
1971  }
1972 
1973  if (NULL != lpResult->pWaveFunctions)
1974  {
1975  for (i = 0; i < lpResult->nEigenValueCount -1; i++)
1976  lpResult->pWaveFunctions[i].Finalize();
1977 
1978  delete[] lpResult->pWaveFunctions;
1979  lpResult->pWaveFunctions = NULL;
1980  lpResult->pWaveFunctions = pWaveFunctions;
1981  }
1982 
1983  if (NULL != lpResult->pEigenVectorsForAMatrix)
1984  {
1985  for (i = 0; i < lpResult->nEigenValueCount-1; i++)
1986  lpResult->pEigenVectorsForAMatrix[i].Finalize();
1987 
1988  delete[] lpResult->pEigenVectorsForAMatrix;
1989  lpResult->pEigenVectorsForAMatrix = NULL;
1990  lpResult->pEigenVectorsForAMatrix = pEigenVectorsForAMatrix;
1991 
1992  if (bInsertFirst)
1993  nNewIndex = 0;
1994  else
1995  nNewIndex = lpResult->nEigenValueCount-1;
1996 
1997  lpResult->pEigenVectorsForAMatrix[nNewIndex].SetSize(pEigenVector->GetSize());
1998  lpResult->pEigenVectorsForAMatrix[nNewIndex] = *pEigenVector;
1999  }
2000 }
2001 
2010 void CKNLanczosMethod::MergeDegeneratedEigenvalues(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, unsigned int nFindingDegeneratedEVCount, CKNMatrixOperation::CKNCSR *pA, CKNMatrixOperation::CKNCSR *pLocalBlock, CKNMatrixOperation::CKNCSR *pLeftBlock, CKNMatrixOperation::CKNCSR *pRightBlock)
2011 {
2012  unsigned int i, j, *pTargetDeflationGroup = NULL, k;
2013  unsigned int *pTargetDeflationEV = NULL, *pEVFindIteration = NULL;
2014  int *pEigenValueCount = NULL;
2015  double *pEVTotal = NULL;
2016  int nEVTotalCount = 0;
2017  MPI_Request req[2];
2018  double Command[COMMAND_SIZE];
2019  bool bKeepWait = true;
2020  unsigned int nTargetGroup;
2021  CKNComplex complexResult;
2022 
2023 
2025  if( CKNMPIManager::IsLanczosComputeRoot() )
2026  {
2027  pEigenValueCount = CKNMPIManager::GetEigenvalueCountFromDeflationGroup(nFindingDegeneratedEVCount, lpResult->nEigenValueCount);
2028 
2029  if( CKNMPIManager::IsDeflationRoot() )
2030  {
2031  int nDeflationStartIndex, nDeflationIndex, nTargetIndex;
2032 
2033  for( i = 0; i < nFindingDegeneratedEVCount ; ++ i)
2034  nEVTotalCount += pEigenValueCount[i];
2035 
2036  pEVTotal = (double*)malloc(sizeof(double)*nEVTotalCount);
2037  pTargetDeflationGroup = (unsigned int*)malloc(sizeof(unsigned int)*nEVTotalCount);
2038  pTargetDeflationEV = (unsigned int*)malloc(sizeof(unsigned int)*nEVTotalCount);
2039  pEVFindIteration = (unsigned int*)malloc(sizeof(unsigned int)*nEVTotalCount);
2040 
2041  nDeflationStartIndex = 0;
2042  nDeflationIndex = 0;
2043 
2044  for( i = 0; i < nFindingDegeneratedEVCount ; ++ i)
2045  {
2046  nTargetIndex = 0;
2047  for( j = nDeflationStartIndex ; j < nDeflationStartIndex + pEigenValueCount[i] ; ++j )
2048  {
2049  pTargetDeflationGroup[j] = nDeflationIndex;
2050  pTargetDeflationEV[j] = nTargetIndex++;
2051  }
2052 
2053  nDeflationStartIndex += pEigenValueCount[i];
2054  nDeflationIndex++;
2055  }
2056  }
2057 
2058  CKNMPIManager::GatherEVFromDeflationGroup(nFindingDegeneratedEVCount, pEVTotal, pEigenValueCount, lpResult->pEigenValues, lpResult->nEigenValueCount);
2059  CKNMPIManager::GatherEVIterationFromDeflationGroup(nFindingDegeneratedEVCount, (int*)pEVFindIteration, pEigenValueCount, (int*)lpResult->pEigenValueFoundIteration, lpResult->nEigenValueCount);
2060  }
2061 
2062  for( i = 0 ; i < lpResult->nEigenValueCount ; ++ i )
2063  FREE_MEM(lpResult->pEigenVectors[i]);
2064  FREE_MEM(lpResult->pEigenVectors);
2065 
2066 
2067  if( CKNMPIManager::IsDeflationRoot() )
2068  {
2069  if( CKNMPIManager::IsLanczosComputeRoot() )
2070  {
2071  int nStartIndex = lpResult->nEigenValueCount;
2072 
2073  lpResult->nDegeneratedEigenValueCount = 0;
2074  lpResult->pDegeneratedEigenValues = (double*)malloc(sizeof(double)*(nEVTotalCount-lpResult->nEigenValueCount));
2075  for( i = pEigenValueCount[0] ; i < nEVTotalCount ; ++i )
2076  {
2077  bool bNewEigenValue = true;
2078 
2080  for( j = 0 ; j < nStartIndex ; ++j )
2081  {
2082  if( CKNMatrixOperation::IsSameA(pEVTotal[i], lpResult->pEigenValues[j], TOLERANE_M_10_9))
2083  {
2084  bool bDoOrthgonal = false;
2085  CKNMatrixOperation::CKNVector vectorFromDeflation;
2086  std::vector<unsigned int> vectorOrthgonalTarget;
2087 
2088  vectorOrthgonalTarget.push_back(j);
2089  bNewEigenValue = false;
2090 
2091  Command[0] = CHECK_EV_ORTH;
2092  Command[1] = (double)pTargetDeflationGroup[i];
2093  Command[2] = (double)pTargetDeflationEV[i];
2094 
2095  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetDeflationComm());
2096  Command[2] = j;
2097  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2098 
2099  vectorFromDeflation.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
2100  CKNMPIManager::ReceiveVectorSync(pTargetDeflationGroup[i], &vectorFromDeflation, vectorFromDeflation.GetSize(), &req[0], CKNMPIManager::GetDeflationComm());
2101 
2102  CKNMatrixOperation::VVDot(&lpResult->pEigenVectorsForAMatrix[j], &vectorFromDeflation, &complexResult);
2103 
2104  bDoOrthgonal = CKNMatrixOperation::IsSameA(fabs(complexResult.GetRealNumber()), 1, GENERAL_TOLERANCE);
2105  bDoOrthgonal = !bDoOrthgonal;
2106 
2107  if( bDoOrthgonal )
2108  {
2109  for( k = nStartIndex ; k < lpResult->nEigenValueCount ; ++k )
2110  {
2111  if( CKNMatrixOperation::IsSameA(pEVTotal[i], lpResult->pEigenValues[k], TOLERANE_M_10_9))
2112  {
2113  Command[0] = CHECK_EV_ORTH_SIMPLE;
2114  Command[1] = (double)k;
2115  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2116 
2117  CKNMatrixOperation::VVDot(&lpResult->pEigenVectorsForAMatrix[k], &vectorFromDeflation, &complexResult);
2118  if( CKNMatrixOperation::IsSameA(fabs(complexResult.GetRealNumber()), 1, GENERAL_TOLERANCE) )
2119  {
2120  bDoOrthgonal = false;
2121  break;
2122  }
2123  else
2124  vectorOrthgonalTarget.push_back(k);
2125  }
2126  }
2127  }
2128 
2129  Command[0] = DO_ORTHGONAL;
2130  Command[2] = j;
2131  Command[1] = bDoOrthgonal ? 1 : 0;
2132 
2133  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2134  if( bDoOrthgonal )
2135  {
2136  int nOrthogonalTarget;
2137  int *pOrthgonalTarget = NULL;
2138  double fEigenValue;
2139 
2140  nOrthogonalTarget = vectorOrthgonalTarget.size();
2141  CKNMPIManager::BroadcastInt(&nOrthogonalTarget, 1);
2142 
2143  pOrthgonalTarget = (int*)malloc(sizeof(int)*nOrthogonalTarget);
2144  for( k = 0 ; k < nOrthogonalTarget ; ++ k )
2145  pOrthgonalTarget[k] = vectorOrthgonalTarget[k];
2146  CKNMPIManager::BroadcastInt(pOrthgonalTarget, nOrthogonalTarget);
2147 
2148  for( k = 0 ; k < nOrthogonalTarget ; ++ k )
2149  CKNMatrixOperation::Gram_schmidt(&lpResult->pEigenVectorsForAMatrix[pOrthgonalTarget[k]], &vectorFromDeflation);
2150 
2151  CKNMatrixOperation::VVDot(&lpResult->pEigenVectorsForAMatrix[j], &vectorFromDeflation, &complexResult);
2152  FREE_MEM(pOrthgonalTarget);
2153 
2154  if( complexResult == 0.0 )
2155  {
2156  CKNMatrixOperation::CKNVector vectorTemp, vectorTemp2;
2157  vectorTemp.SetSize(vectorFromDeflation.GetSize());
2158  vectorTemp2.SetSize(vectorFromDeflation.GetSize());
2159 
2160  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MVMUL);
2161  CKNMatrixOperation::MVMul(pA, &vectorFromDeflation, &vectorTemp);
2162  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MVMUL);
2163 
2164  vectorTemp2 = vectorFromDeflation;
2165 
2166  if( CKNMPIManager::IsRootRank() )
2167  fEigenValue = lpResult->pEigenValues[j];
2168  CKNMPIManager::BroadcastDouble(&fEigenValue, 1);
2169 
2170  vectorTemp2.ScalarMultiple(fEigenValue);
2171  vectorTemp2.MinusVector(&vectorTemp);
2172 
2173  double fNorm = vectorTemp2.GetNorm(true);
2174 
2175  if( CKNMatrixOperation::IsSame(fNorm, 0.0, GENERAL_TOLERANCE))
2176  {
2177  Command[0] = SEND_BACK_EV;
2178  Command[1] = (double)pTargetDeflationGroup[i];
2179  Command[2] = lpResult->nDegeneratedEigenValueCount+1;
2180  lpResult->pDegeneratedEigenValues[lpResult->nDegeneratedEigenValueCount++] = complexResult.GetRealNumber();
2181  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetDeflationComm());
2182  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2183  CKNMPIManager::SendVectorSync(pTargetDeflationGroup[i], &vectorFromDeflation, vectorFromDeflation.GetSize(), &req[1], CKNMPIManager::GetDeflationComm());
2184  AppendEigenValue(lpResult, pEVTotal[i]);
2185  AppendEigenVector(lpResult, &vectorFromDeflation);
2186  break;
2187  }
2188  else
2189  {
2190  Command[0] = NOT_SEND_BACK_EV;
2191  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2192  }
2193  }
2194  else
2195  {
2196  Command[0] = NOT_SEND_BACK_EV;
2197  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2198  }
2199  }
2200  vectorFromDeflation.Finalize();
2201  }
2202  }
2203 
2204  if(bNewEigenValue)
2205  {
2206  CKNMatrixOperation::CKNVector vectorFromDeflation;
2207 
2208  Command[0] = SEND_EV_TO_MASTER;
2209  Command[1] = (double)pTargetDeflationGroup[i];
2210  Command[2] = (double)pTargetDeflationEV[i];
2211 
2212  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetDeflationComm());
2213  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2214 
2215  vectorFromDeflation.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
2216  CKNMPIManager::ReceiveVectorSync(pTargetDeflationGroup[i], &vectorFromDeflation, vectorFromDeflation.GetSize(), &req[0], CKNMPIManager::GetDeflationComm());
2217 
2218  AppendEigenValue(lpResult, pEVTotal[i], pEVFindIteration[i], true);
2219  AppendEigenVector(lpResult, &vectorFromDeflation, true);
2220 
2221  nStartIndex++;
2222  }
2223  }
2224 
2225  lpResult->pDegeneratedEigenValues = (double*)realloc(lpResult->pDegeneratedEigenValues, sizeof(double)*lpResult->nDegeneratedEigenValueCount);
2226 
2227  Command[0] = EXIT_MPI_WAIT;
2228  Command[1] = -1;
2229  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetDeflationComm());
2230  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2231  }
2232  else
2233  {
2234  while(bKeepWait)
2235  {
2236  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2237 
2238  switch((int)Command[0])
2239  {
2240  case EXIT_MPI_WAIT:
2241  bKeepWait = false;
2242  break;
2243  case SEND_EV_TO_MASTER:
2244  {
2245  CKNMatrixOperation::CKNVector vectorFromDeflation;
2246  vectorFromDeflation.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
2247 
2248  nTargetGroup = (unsigned int)Command[1];
2249  nTargetGroup = CKNMPIManager::GetTotalNodeCount() * nTargetGroup + CKNMPIManager::GetCurrentRank();
2250 
2251  CKNMPIManager::ReceiveVectorSync(nTargetGroup, &vectorFromDeflation, vectorFromDeflation.GetSize(), &req[0], MPI_COMM_WORLD);
2252 
2253  lpResult->nEigenValueCount++;
2254  AppendEigenVector(lpResult, &vectorFromDeflation, true);
2255  }
2256  break;
2257  case CHECK_EV_ORTH:
2258  {
2259  CKNMatrixOperation::CKNVector vectorFromDeflation;
2260  vectorFromDeflation.SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
2261 
2262  nTargetGroup = (unsigned int)Command[1];
2263  nTargetGroup = CKNMPIManager::GetTotalNodeCount() * nTargetGroup + CKNMPIManager::GetCurrentRank();
2264 
2265  CKNMPIManager::ReceiveVectorSync(nTargetGroup, &vectorFromDeflation, vectorFromDeflation.GetSize(), &req[0], MPI_COMM_WORLD);
2266  CKNMatrixOperation::VVDot(&lpResult->pEigenVectorsForAMatrix[(int)Command[2]], &vectorFromDeflation, &complexResult);
2267 
2268  while(1)
2269  {
2270  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2271  if( DO_ORTHGONAL == (int)Command[0] )
2272  break;
2273  else if( CHECK_EV_ORTH_SIMPLE == (int)Command[0] )
2274  CKNMatrixOperation::VVDot(&lpResult->pEigenVectorsForAMatrix[(int)Command[1]], &vectorFromDeflation, &complexResult);
2275  }
2276 
2277  if( DO_ORTHGONAL == (int)Command[0] && 1 == (int)Command[1])
2278  {
2279 
2280  int nOrthogonalTarget;
2281  int *pOrthgonalTarget = NULL;
2282  double fEigenValue;
2283 
2284  CKNMPIManager::BroadcastInt(&nOrthogonalTarget, 1);
2285 
2286  pOrthgonalTarget = (int*)malloc(sizeof(int)*nOrthogonalTarget);
2287  CKNMPIManager::BroadcastInt(pOrthgonalTarget, nOrthogonalTarget);
2288 
2289 
2290  for( k = 0 ; k < nOrthogonalTarget ; ++ k )
2291  CKNMatrixOperation::Gram_schmidt(&lpResult->pEigenVectorsForAMatrix[pOrthgonalTarget[k]], &vectorFromDeflation);
2292 
2293  CKNMatrixOperation::VVDot(&lpResult->pEigenVectorsForAMatrix[(int)Command[2]], &vectorFromDeflation, &complexResult);
2294  FREE_MEM(pOrthgonalTarget);
2295 
2296  if( complexResult == 0.0 )
2297  {
2298  CKNMatrixOperation::CKNVector vectorTemp, vectorTemp2;
2299  vectorTemp.SetSize(vectorFromDeflation.GetSize());
2300  vectorTemp2.SetSize(vectorFromDeflation.GetSize());
2301 
2302  CKNTimeMeasurement::MeasurementStart(CKNTimeMeasurement::MVMUL);
2303  CKNMatrixOperation::MVMul(pA, &vectorFromDeflation, &vectorTemp);
2304  CKNTimeMeasurement::MeasurementEnd(CKNTimeMeasurement::MVMUL);
2305 
2306  vectorTemp2 = vectorFromDeflation;
2307  CKNMPIManager::BroadcastDouble(&fEigenValue, 1);
2308  vectorTemp2.ScalarMultiple(fEigenValue);
2309  vectorTemp2.MinusVector(&vectorTemp);
2310  double fNorm = vectorTemp2.GetNorm(true);
2311  }
2312 
2313  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2314  if( SEND_BACK_EV == (int)Command[0] )
2315  {
2316  CKNMPIManager::SendVectorSync(nTargetGroup, &vectorFromDeflation, vectorFromDeflation.GetSize(), &req[1], MPI_COMM_WORLD);
2317  lpResult->nEigenValueCount++;
2318  lpResult->nDegeneratedEigenValueCount++;
2319  AppendEigenVector(lpResult, &vectorFromDeflation);
2320  }
2321  }
2322 
2323  vectorFromDeflation.Finalize();
2324  }
2325  break;
2326  }
2327  }
2328  }
2329  }
2330  else
2331  {
2332  int nSentEVIndex = -1;
2333 
2334  lpResult->pDegeneratedIndex = (int*)malloc(sizeof(int)*lpResult->nEigenValueCount);
2335  for( i = 0 ; i < lpResult->nEigenValueCount ; ++i)
2336  lpResult->pDegeneratedIndex[i] = -1;
2337 
2338  if( CKNMPIManager::IsLanczosComputeRoot() )
2339  {
2341  while(bKeepWait)
2342  {
2343  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetDeflationComm());
2344  switch((int)Command[0])
2345  {
2346  case EXIT_MPI_WAIT:
2347  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2348  bKeepWait = false;
2349  break;
2350  case CHECK_EV_ORTH:
2351  case SEND_EV_TO_MASTER:
2352  nTargetGroup = (unsigned int)Command[1];
2353  if( nTargetGroup == CKNMPIManager::GetLanczosGroupIndex())
2354  {
2355  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2356  nSentEVIndex = (int)Command[2];
2357  CKNMPIManager::SendVectorSync(0, &lpResult->pEigenVectorsForAMatrix[nSentEVIndex], CKNMPIManager::GetCurrentLoadBalanceCount(), &req[0], CKNMPIManager::GetDeflationComm());
2358  }
2359  break;
2360  case SEND_BACK_EV:
2361  nTargetGroup = (unsigned int)Command[1];
2362  if( nTargetGroup == CKNMPIManager::GetLanczosGroupIndex())
2363  {
2364  lpResult->pDegeneratedIndex[nSentEVIndex] = (unsigned int)Command[2];
2365  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2366  CKNMPIManager::ReceiveVectorSync(0, &lpResult->pEigenVectorsForAMatrix[nSentEVIndex], CKNMPIManager::GetCurrentLoadBalanceCount(), &req[1], CKNMPIManager::GetDeflationComm());
2367  }
2368  break;
2369  }
2370  }
2371  }
2372  else
2373  {
2374  while(bKeepWait)
2375  {
2376  CKNMPIManager::ExchangeCommand(Command, CKNMPIManager::GetLanczosComputComm());
2377  switch((int)Command[0])
2378  {
2379  case EXIT_MPI_WAIT:
2380  bKeepWait = false;
2381  break;
2382  case SEND_EV_TO_MASTER:
2383  case CHECK_EV_ORTH:
2384  {
2385  int nLanczosGroupSize = CKNMPIManager::GetTotalNodeCount();
2386  nSentEVIndex = (int)Command[2];
2387  CKNMPIManager::SendVectorSync(CKNMPIManager::GetCurrentRank(), &lpResult->pEigenVectorsForAMatrix[nSentEVIndex], CKNMPIManager::GetCurrentLoadBalanceCount(), &req[0], MPI_COMM_WORLD);
2388  }
2389  break;
2390  case SEND_BACK_EV:
2391  lpResult->pDegeneratedIndex[nSentEVIndex] = (unsigned int)Command[2];
2392  CKNMPIManager::ReceiveVectorSync(CKNMPIManager::GetCurrentRank(), &lpResult->pEigenVectorsForAMatrix[nSentEVIndex], CKNMPIManager::GetCurrentLoadBalanceCount(), &req[1], MPI_COMM_WORLD);
2393  break;
2394  }
2395  }
2396  }
2397  }
2398 
2399  CKNMPIManager::BarrierAllComm();
2400 
2401  FREE_MEM(pEVTotal);
2402  FREE_MEM(pEigenValueCount);
2403  FREE_MEM(pEVFindIteration);
2404  FREE_MEM(pTargetDeflationGroup);
2405 }
2406 
2413 bool CKNLanczosMethod::CheckingCalculationCondition(bool bCalcuEigenValue, bool bCalcuWaveFunction, unsigned int nDeflationGroup)
2414 {
2415  bool bRtn = false;
2416 
2417  if( true == bCalcuWaveFunction )
2418  if( false == bCalcuEigenValue )
2419  return bRtn;
2420 
2421  if( nDeflationGroup > 1 && false == bCalcuEigenValue )
2422  return bRtn;
2423 
2424  bRtn = true;
2425  return bRtn;
2426 }
2427 
2433 int CKNLanczosMethod::ResultCompare(const void *pA, const void *pB)
2434 {
2435  CKNLanczosMethod::LPRESULT_SORT_DATA lpA = (CKNLanczosMethod::LPRESULT_SORT_DATA)pA;
2436  CKNLanczosMethod::LPRESULT_SORT_DATA lpB = (CKNLanczosMethod::LPRESULT_SORT_DATA)pB;
2437 
2438 
2439  if (lpA->fEigenValue > lpB->fEigenValue)
2440  return 1;
2441  else if(lpA->fEigenValue < lpB->fEigenValue)
2442  return -1;
2443  else
2444  {
2445  if( (unsigned int)lpA->nEigenValueFoundIteration > (unsigned int)lpB->nEigenValueFoundIteration )
2446  return 1;
2447  else
2448  return -1;
2449  }
2450 
2451  return -1;
2452 }
2453 
2457 void CKNLanczosMethod::SortSolution(LPEIGENVALUE_RESULT lpResult)
2458 {
2459 
2460  LPRESULT_SORT_DATA lpData = NULL;
2461  int *pOrder = NULL;
2462  unsigned int i;
2463  CKNMatrixOperation::CKNVector *pVectorEV = NULL, *pVectorWF = NULL;
2464 
2465  if( lpResult->nEigenValueCount <= 1 || false == CKNMPIManager::IsRootRank() )
2466  return;
2467 
2468  pOrder = (int*)malloc(sizeof(int)*lpResult->nEigenValueCount*2);
2469  if( NULL != lpResult->pEigenVectorsForAMatrix )
2470  pVectorEV = new CKNMatrixOperation::CKNVector[lpResult->nEigenValueCount];
2471  if( NULL != lpResult->pWaveFunctions )
2472  pVectorWF = new CKNMatrixOperation::CKNVector[lpResult->nEigenValueCount];
2473 
2474  if( CKNMPIManager::IsRootRank() )
2475  {
2476  lpData = (LPRESULT_SORT_DATA)malloc(sizeof(RESULT_SORT_DATA)*lpResult->nEigenValueCount);
2477 
2478  for( i = 0 ; i < lpResult->nEigenValueCount ; ++ i )
2479  {
2480  lpData[i].fEigenValue = lpResult->pEigenValues[i];
2481  lpData[i].nEigenValueFoundIteration = lpResult->pEigenValueFoundIteration[i];
2482  lpData[i].nOriginalIndex = i;
2483  }
2484 
2485  qsort(lpData, lpResult->nEigenValueCount, sizeof(RESULT_SORT_DATA), CKNLanczosMethod::ResultCompare);
2486 
2487  for( i = 0 ; i < lpResult->nEigenValueCount ; ++ i )
2488  {
2489  pOrder[i*2] = i;
2490  pOrder[i*2+1] = lpData[i].nOriginalIndex;
2491  lpResult->pEigenValues[i] = lpData[i].fEigenValue;
2492  lpResult->pEigenValueFoundIteration[i] = lpData[i].nEigenValueFoundIteration;
2493  }
2494  }
2495 
2496  /*if( CKNMPIManager::IsMultiLevelMPI() )
2497  CKNMPIManager::BroadcastInt(pOrder, lpResult->nEigenValueCount*2);*/
2498 
2499  for( i = 0; i < lpResult->nEigenValueCount ; ++ i)
2500  {
2501  if( NULL != lpResult->pEigenVectorsForAMatrix )
2502  {
2503  pVectorEV[i].SetSize(CKNMPIManager::GetCurrentLoadBalanceCount());
2504  pVectorEV[i] = lpResult->pEigenVectorsForAMatrix[i];
2505  }
2506  if( NULL != lpResult->pWaveFunctions )
2507  {
2508  pVectorWF[i].SetSize(CKNMPIManager::GetCurrentLoadBalanceCount()/10);
2509  pVectorWF[i] = lpResult->pWaveFunctions[i];
2510  }
2511  }
2512 
2513  for( i = 0; i < lpResult->nEigenValueCount ; ++ i)
2514  {
2515  if( NULL != lpResult->pEigenVectorsForAMatrix )
2516  lpResult->pEigenVectorsForAMatrix[pOrder[i*2]] = pVectorEV[pOrder[i*2+1]];
2517  if( NULL != lpResult->pWaveFunctions )
2518  lpResult->pWaveFunctions[pOrder[i*2]] = pVectorWF[pOrder[i*2+1]];
2519  }
2520 
2521  for( i = 0; i < lpResult->nEigenValueCount ; ++ i)
2522  {
2523  if( NULL != lpResult->pEigenVectorsForAMatrix )
2524  pVectorEV[i].Finalize();
2525  if( NULL != lpResult->pWaveFunctions )
2526  pVectorWF[i].Finalize();
2527  }
2528  if( NULL != lpResult->pEigenVectorsForAMatrix )
2529  delete[] pVectorEV;
2530  if( NULL != lpResult->pWaveFunctions )
2531  delete[] pVectorWF;
2532 
2533  FREE_MEM(lpData);
2534  FREE_MEM(pOrder);
2535 }
CKNMPIManager::BroadcastBool
static void BroadcastBool(bool *boolValue, int nRootRank=0)
Broadcst boolean value.
Definition: KNMPIManager.cpp:524
CKNMPIManager::GatherEVFromDeflationGroup
static void GatherEVFromDeflationGroup(int nSourceCount, double *pReceiveBuffer, int *pSourceCount, double *pSendBuffer, int nSendCount)
Definition: KNMPIManager.h:76
CKNMatrixOperation::CKNVector::SetSize
void SetSize(unsigned int nSize)
Set Vector elements size.
Definition: KNMatrixOperation.cpp:40
DO_ORTHGONAL
#define DO_ORTHGONAL
Definition: CKNGlobal.h:96
CKNLanczosMethod::EIGENVALUE_RESULT::pWaveFunctions
CKNMatrixOperation::CKNVector * pWaveFunctions
Definition: KNLanczosMethod.h:40
CKNLanczosMethod::FinalLanczosVector
void FinalLanczosVector()
Deallocating lanczos vectors.
Definition: KNLanczosMethod.cpp:595
CKNComplex::GetImaginaryNumber
double GetImaginaryNumber() const
Get imaginary part.
Definition: KNComplex.h:27
CKNMPIManager::BroadcastInt
static void BroadcastInt(int *pValue, unsigned int nSize, int nRootRank=0, MPI_Comm comm=MPI_COMM_NULL)
Broadcst boolean value.
Definition: KNMPIManager.cpp:560
CKNMatrixOperation::CKNVector::ScalarDivision
void ScalarDivision(CKNComplex Scalar)
Scalar division operation.
Definition: KNMatrixOperation.cpp:164
CKNIPCCUtility::ShowMsg
static void ShowMsg(char *pszBuffer)
Show message.
Definition: KNIPCCUtility.cpp:35
CKNMatrixOperation::CKNVector::ScalarMultiple
void ScalarMultiple(CKNComplex Scalar)
Scalar multiple operation.
Definition: KNMatrixOperation.cpp:129
CKNMatrixOperation::CKNVector::Normalize
void Normalize(bool bMPI=false)
Normalize vector with norm.
Definition: KNMatrixOperation.cpp:216
CKNMPIManager::Barrier
static void Barrier()
Definition: KNMPIManager.h:70
CKNLanczosMethod::DoSelectiveReorthogonalization
void DoSelectiveReorthogonalization(unsigned int nIterationCount)
Do selective reorthogonalization.
Definition: KNLanczosMethod.cpp:719
CKNMatrixOperation::CKNVector::ScalarMultiThanMinusVector
void ScalarMultiThanMinusVector(double fScalar, CKNVector *vector)
Do minus operation after scalar multiple to operand between vectors.
Definition: KNMatrixOperation.cpp:264
CKNMPIManager::GetMPIComm
static MPI_Comm GetMPIComm()
Definition: KNMPIManager.h:67
X_largest
unsigned int X_largest
Definition: KNLanczosMethod.cpp:27
CKNLanczosMethod::EIGENVALUE_RESULT::pEigenValues
double * pEigenValues
Definition: KNLanczosMethod.h:31
CKNLanczosMethod::MergeDegeneratedEigenvalues
static void MergeDegeneratedEigenvalues(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, unsigned int nFindingDegeneratedEVCount, CKNMatrixOperation::CKNCSR *pA, CKNMatrixOperation::CKNCSR *pLocalBlock, CKNMatrixOperation::CKNCSR *pLeftBlock, CKNMatrixOperation::CKNCSR *pRightBlock)
Merging eigenvalue into mater group.
Definition: KNLanczosMethod.cpp:2010
CKNMPIManager::GetEigenvalueCountFromDeflationGroup
static int * GetEigenvalueCountFromDeflationGroup(int nDeflationGroupCount, int nLocalEVCount)
Checking is root rank of Deflation computation.
Definition: KNMPIManager.cpp:802
CKNMPIManager::BroadcastLanczosResult
static void BroadcastLanczosResult(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, int nIterationCount)
Broadcast Lanczos result.
Definition: KNMPIManager.cpp:578
CKNLanczosMethod::DoResidualCheck
void DoResidualCheck(CKNMatrixOperation::CKNCSR *pAMatrix, LPEIGENVALUE_RESULT lpResult)
Residual checking.
Definition: KNLanczosMethod.cpp:1565
KNLanczosMethod.h
Show message and debugging variable.
CKNMPIManager::BroadcastDouble
static void BroadcastDouble(double *pValue, unsigned int nSize, int nRootRank=0, MPI_Comm comm=MPI_COMM_NULL)
Broadcst boolean value.
Definition: KNMPIManager.cpp:540
CKNComplex::GetRealNumber
double GetRealNumber() const
Get real part.
Definition: KNComplex.h:26
CKNLanczosMethod::EIGENVALUE_RESULT::pEigenVectorsForAMatrix
CKNMatrixOperation::CKNVector * pEigenVectorsForAMatrix
Definition: KNLanczosMethod.h:35
CKNLanczosMethod::InitLanczosIterationVariables
void InitLanczosIterationVariables(CKNComplex **pAlpha, double **pAlphaReal, double **pBeta, double **pWj, double **pWjm1, double **pWjp1, CKNMatrixOperation::CKNVector **pW)
Init omega, alpha, beta array.
Definition: KNLanczosMethod.cpp:617
CKNMatrixOperation::CKNCSR::GetColumnCount
unsigned int GetColumnCount()
Getting row size of matrix.
Definition: KNMatrixOperation.h:173
CKNMPIManager::GetLanczosComputComm
static MPI_Comm GetLanczosComputComm()
Definition: KNMPIManager.h:79
CKNMatrixOperation::CKNCSR
Data and operation representation of CSR(Compressed Sparse Row)
Definition: KNMatrixOperation.h:164
CKNMPIManager::GetCurrentLoadBalanceCount
static int GetCurrentLoadBalanceCount()
Get Current node's rank load balancing number.
Definition: KNMPIManager.cpp:684
CKNLanczosMethod::LPRESULT_SORT_DATA
struct CKNLanczosMethod::RESULT_SORT_DATA * LPRESULT_SORT_DATA
CKNLanczosMethod::InitLanczosVector
void InitLanczosVector()
Init lanczos vectors.
Definition: KNLanczosMethod.cpp:576
CKNTimeMeasurement::MeasurementEnd
static void MeasurementEnd(MEASUREMENT_INDEX index)
Measurement end for part.
Definition: KNTimeMeasurement.cpp:87
CKNLanczosMethod::EIGENVALUE_RESULT::nEigenVectorSize
unsigned int nEigenVectorSize
Definition: KNLanczosMethod.h:34
CKNLanczosMethod::ResultCompare
static int ResultCompare(const void *pA, const void *pB)
Comparing computing result function for quick sorting.
Definition: KNLanczosMethod.cpp:2433
CKNTimeMeasurement::GetTakeTime
static double GetTakeTime(MEASUREMENT_INDEX index)
Get taken time for part.
Definition: KNTimeMeasurement.cpp:61
CKNMatrixOperation::MVMulEx_Optimal
static void MVMulEx_Optimal(CKNCSR *pAMatrix, CKNVector *pVector, CKNVector *pResult, unsigned int, unsigned int, CKNVector *, int)
Matrix and vector multiple operation for 1 layer exchanging communication.
Definition: KNMatrixOperation.cpp:2452
CKNMPIManager::GetLanczosGroupIndex
static unsigned int GetLanczosGroupIndex()
Definition: KNMPIManager.h:83
CKNMPIManager::GatherEVIterationFromDeflationGroup
static void GatherEVIterationFromDeflationGroup(int nSourceCount, int *pReceiveBuffer, int *pSourceCount, int *pSendBuffer, int nSendCount)
Gather eigenvalue from deflation group.
Definition: KNMPIManager.h:77
CKNLanczosMethod::InitializeTemporaryArrayAndVector
bool InitializeTemporaryArrayAndVector(int nIterationCount)
Initialize temporary eigenvalue arrays and vectors.
Definition: KNLanczosMethod.cpp:76
CKNLanczosMethod::EIGENVALUE_RESULT::pDegeneratedIndex
int * pDegeneratedIndex
Definition: KNLanczosMethod.h:42
CKNMPIManager::GetTotalNodeCount
static int GetTotalNodeCount()
Definition: KNMPIManager.h:44
CKNLanczosMethod::BuildWaveFunction
void BuildWaveFunction(LPEIGENVALUE_RESULT lpResult)
Building wavefunction.
Definition: KNLanczosMethod.cpp:1530
CKNMPIManager::ReceiveVectorSync
static void ReceiveVectorSync(int nSourceRank, CKNMatrixOperation::CKNVector *pVector, int nSize, MPI_Request *req, MPI_Comm commWorld=MPI_COMM_NULL)
Receiving Vector with sync.
Definition: KNMPIManager.cpp:910
SHOW_SIMPLE_MSG
#define SHOW_SIMPLE_MSG(message)
Definition: CKNGlobal.h:41
CKNLanczosMethod::FinalizeLanczosInterationVariable
void FinalizeLanczosInterationVariable(CKNComplex *pAlpha, double *pAlphaReal, double *pBeta, double *pWj, double *pWjm1, double *pWjp1, CKNMatrixOperation::CKNVector *pW)
Deallocating omega, alpha, beta.
Definition: KNLanczosMethod.cpp:645
XeonPhi_header.h
CKNLanczosMethod::RESULT_SORT_DATA::fEigenValue
double fEigenValue
Definition: KNLanczosMethod.h:50
GENERAL_TOLERANCE
#define GENERAL_TOLERANCE
General tolerance definition.
Definition: CKNGlobal.h:48
CKNLanczosMethod::CKNLanczosMethod
CKNLanczosMethod()
Constructor.
Definition: KNLanczosMethod.cpp:32
CKNMatrixOperation::CKNVector::MinusVector
void MinusVector(CKNVector *vector)
Do minus operation between vectors.
Definition: KNMatrixOperation.cpp:322
CKNLanczosMethod::RecalcuWaveFunction
static void RecalcuWaveFunction(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult)
Recalculating wavefunction after merging degenerated eigenvalues.
Definition: KNLanczosMethod.cpp:1873
CKNLanczosMethod::CheckAndDoSelectiveReorthogonalization
bool CheckAndDoSelectiveReorthogonalization(int nIterationCount, double *pAlpha, double *pBeta, double *pWj, double *pWjm1, double *pWjp1, double fANorm)
Check current state need selective reorthogonalization and do it.
Definition: KNLanczosMethod.cpp:708
TOLERANE_M_10_9
#define TOLERANE_M_10_9
10^-9
Definition: CKNGlobal.h:49
CKNLanczosMethod::AppendEigenVector
static void AppendEigenVector(LPEIGENVALUE_RESULT lpResult, CKNMatrixOperation::CKNVector *pEigenVector, bool bInsertFirst=false)
Appending eigenvector into master group if degenerated eigenvalue is finded.
Definition: KNLanczosMethod.cpp:1941
CKNLanczosMethod::IntegrateEigenvalues
void IntegrateEigenvalues(int nIterationCount, LPEIGENVALUE_RESULT lpResult, unsigned int nCalculatedEigenValueCount, double *pCalcuResult_Value, double *pCalcuResult_Vector)
Integrating computing solution during Lanczos method operation.
Definition: KNLanczosMethod.cpp:1051
CKNMatrixOperation::IsSameA
static bool IsSameA(double operand1, double operand2, double tol)
Compare two double variable.
Definition: KNMatrixOperation.cpp:2802
CKNLanczosMethod::RESULT_SORT_DATA
Structure for eigenvalue result sorting.
Definition: KNLanczosMethod.h:49
CKNLanczosMethod::BuildTMatrix
double * BuildTMatrix(unsigned int nOrder, double *pAlpha, double *pBeta)
Building T matrix for solving eigenvalue.
Definition: KNLanczosMethod.cpp:1367
CKNLanczosMethod::StopIteration
static void StopIteration()
Stop lanczos iteration on going state.
Definition: KNLanczosMethod.cpp:1522
DEGENERATED_INDEX
#define DEGENERATED_INDEX
Definition: CKNGlobal.h:56
CKNLanczosMethod::SaveLanczosResult
static void SaveLanczosResult(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, bool bCalcuEigenvalue, bool bWaveFunction, double *pKValue, int nRepeatCount)
Saving Lanczos computation result into file.
Definition: KNLanczosMethod.cpp:1654
EXIT_MPI_WAIT
#define EXIT_MPI_WAIT
Command for Deflation Lanczos.
Definition: CKNGlobal.h:94
CKNGlobal.h
Common definition for Solver.
CKNMPIManager::SendVectorSync
static void SendVectorSync(int nTargetRank, CKNMatrixOperation::CKNVector *pVector, int nSize, MPI_Request *req, MPI_Comm commWorld=MPI_COMM_NULL)
Getting Deflation computing group MPI_Comm.
Definition: KNMPIManager.cpp:890
CKNLanczosMethod::CheckingCalculationCondition
static bool CheckingCalculationCondition(bool bCalcuEigenValue, bool bCalcuWaveFunction, unsigned int nDeflationGroup)
Checking pre conditions for Lanczos method operation.
Definition: KNLanczosMethod.cpp:2413
CKNMatrixOperation::CKNVector::GetAt
CKNComplex GetAt(unsigned int nIndex)
Get element value from specific index.
Definition: KNMatrixOperation.cpp:100
CHECK_EV_ORTH
#define CHECK_EV_ORTH
Definition: CKNGlobal.h:95
CKNMPIManager::IsLanczosComputeRoot
static bool IsLanczosComputeRoot()
Barrier current deflation group.
Definition: KNMPIManager.h:71
CKNLanczosMethod::EIGENVALUE_RESULT
Structure for engienvalue computing.
Definition: KNLanczosMethod.h:30
CKNLanczosMethod::DoEigenValueSolving
bool DoEigenValueSolving(int nIterationCount, double *pAlpha, double *pBeta, double fANorm, LPEIGENVALUE_RESULT lpResult, bool bFinal)
Every user set iteration count calculate eigenvalues.
Definition: KNLanczosMethod.cpp:872
CKNLanczosMethod::ConvergenceCheckingEx
int ConvergenceCheckingEx(int nEigenValueCount, double *pEigenValues, double *pEiegnVectors, bool *pbValidEigenValue, double fANorm, double *pBeta, int nIterationCount)
Checking convergence criteria.
Definition: KNLanczosMethod.cpp:1101
CHECK_EV_ORTH_SIMPLE
#define CHECK_EV_ORTH_SIMPLE
Definition: CKNGlobal.h:100
CKNMatrixOperation::MVMul
static void MVMul(CKNCSR *pAMatrix, CKNVector *pVector, CKNVector *pResult)
Matrix and vector multiple operation.
Definition: KNMatrixOperation.cpp:2134
CKNLanczosMethod::EIGENVALUE_RESULT::nDegeneratedEigenValueCount
unsigned int nDegeneratedEigenValueCount
Definition: KNLanczosMethod.h:37
CKNMatrixOperation::MVMulOptimal
static void MVMulOptimal(CKNCSR *pAMatrix, CKNVector *pVector, CKNVector *pResult)
Matrix and vector multiple operation for multiple call.
Definition: KNMatrixOperation.cpp:2600
CKNLanczosMethod::EIGENVALUE_RESULT::pEigenValueFoundIteration
unsigned int * pEigenValueFoundIteration
Definition: KNLanczosMethod.h:39
KNTimeMeasurement.h
Time measurement class.
CKNLanczosMethod::SpuriousRitzValueChecking
int SpuriousRitzValueChecking(int nEigenValueCount, double *pEigenValues, double *pEigenVectors, double *pNonSpuriousValues, double *pNonSpuriousVectors, double fANorm, int nIterationCount)
Checking spurious values.
Definition: KNLanczosMethod.cpp:1202
SEND_EV_TO_MASTER
#define SEND_EV_TO_MASTER
Definition: CKNGlobal.h:99
CKNMatrixOperation::CKNVector::m_vectValueImaginaryBuffer
double_vector_t m_vectValueImaginaryBuffer
A member variable for saving none zero elements.
Definition: KNMatrixOperation.h:105
CKNLanczosMethod::EIGENVALUE_RESULT::nMaxEigenValueFoundIteration
unsigned int nMaxEigenValueFoundIteration
Definition: KNLanczosMethod.h:38
CKNLanczosMethod::ExtractDoubleVector
void ExtractDoubleVector(unsigned int nVectorsize, double *pTarget, double *pSource, unsigned int nSrcCount, int *pFilter, unsigned int nFilterCount, bool bExclusive)
Extract vectors by condition that described in filter.
Definition: KNLanczosMethod.cpp:1453
inverse_iter
void inverse_iter(double *alpha, double *beta, double *app_evc, int iter, double app_eva)
Definition: KNLanczosMethod.cpp:790
NOT_SEND_BACK_EV
#define NOT_SEND_BACK_EV
Definition: CKNGlobal.h:98
CKNLanczosMethod::EIGENVALUE_RESULT::pDegeneratedEigenValues
double * pDegeneratedEigenValues
Definition: KNLanczosMethod.h:32
CKNMatrixOperation::CKNVector::Finalize
void Finalize()
Free allocated memory for vector elements.
Definition: KNMatrixOperation.cpp:365
CKNComplex::GetAbsolute
double GetAbsolute()
Get Absolute value of complex number.
Definition: KNComplex.cpp:24
CKNLanczosMethod::DoLanczosMethod
LPEIGENVALUE_RESULT DoLanczosMethod(CKNMatrixOperation::CKNCSR *pAMatrix, unsigned int nIterationCount, unsigned int nEigenValueCheckInterval, unsigned int nEigenValueCount, double fEigenvalueMin, double fEignevalueMax, double fConvergenceTolerance, bool bReorthogonalization, bool bCalcuEigVector, bool bWaveFunction, double load_in_MIC, CKNMatrixOperation::CKNCSR *pmylocalblock=NULL, CKNMatrixOperation::CKNCSR *leftlocalblock=NULL, CKNMatrixOperation::CKNCSR *rightlocalblock=NULL)
Doing lanczos method.
Definition: KNLanczosMethod.cpp:130
CKNLanczosMethod::InitVariables
void InitVariables()
Deallocating member variables.
Definition: KNLanczosMethod.cpp:1395
KNIPCCUtility.h
This class includes functions for matrix debugging.
CKNTimeMeasurement::GetTotalTakeTime
static double GetTotalTakeTime()
Definition: KNTimeMeasurement.h:32
CKNLanczosMethod::SortSolution
static void SortSolution(LPEIGENVALUE_RESULT lpResult)
Sorting computing eigenvalue.
Definition: KNLanczosMethod.cpp:2457
CKNLanczosMethod::DistinguishClusterOfEigenvalue
int DistinguishClusterOfEigenvalue(int nEigenValueCount, double *pEigenValues, double *pEigenVectors, double *pNonClustersValues, double *pNonClustersVectors, int nIterationCount)
Distinguish clusters values.
Definition: KNLanczosMethod.cpp:1271
CKNMPIManager::ExchangeCommand
static void ExchangeCommand(double *pfCommand, MPI_Comm comm)
Gather eigenvalue finding iteration number from deflation group.
Definition: KNMPIManager.cpp:878
CKNMatrixOperation::CKNVector::GetNorm
double GetNorm(bool bMPI=false)
Getting norm of vector.
Definition: KNMatrixOperation.cpp:232
CKNLanczosMethod::CalculateEigenVector
void CalculateEigenVector(LPEIGENVALUE_RESULT lpResult, CKNMatrixOperation::CKNVector V, unsigned int nIterationIndex)
Calculate Eigen vector of A Matrix.
Definition: KNLanczosMethod.cpp:673
CKNMPIManager::GetCurrentRank
static int GetCurrentRank()
Definition: KNMPIManager.h:42
CKNLanczosMethod::RESULT_SORT_DATA::nOriginalIndex
unsigned int nOriginalIndex
Definition: KNLanczosMethod.h:52
thomas_alg
void thomas_alg(double **T, double *initguess, double *app_evc, int iter)
Definition: KNLanczosMethod.cpp:725
CKNMPIManager::GetDeflationComm
static MPI_Comm GetDeflationComm()
Getting Lanczos computing group MPI_Comm.
Definition: KNMPIManager.h:80
CKNLanczosMethod::IntegrateEigenvaluesEx
void IntegrateEigenvaluesEx(int nIterationCount, LPEIGENVALUE_RESULT lpResult, unsigned int nCalculatedEigenValueCount, unsigned int nCalculatedEigenValueCountBeforeConvergenceCheck, double *pCalcuResult_Value, double *pCalcuResult_Vector, bool *pbValidEigenValue)
Integrating computing solution during Lanczos method operation.
Definition: KNLanczosMethod.cpp:948
CKNMatrixOperation::CKNVector::m_vectValueRealBuffer
double_vector_t m_vectValueRealBuffer
A member variable for saving none zero elements.
Definition: KNMatrixOperation.h:104
ALLOC_WITH_NULL_INIT
#define ALLOC_WITH_NULL_INIT(pointer, data_type, data_size)
Definition: CKNGlobal.h:35
CKNLanczosMethod::AppendEigenValue
static void AppendEigenValue(LPEIGENVALUE_RESULT lpResult, double fEigenValue, unsigned int nFindIteration=DEGENERATED_INDEX, bool bInsertFirst=false)
Checking is aborting computation flag.
Definition: KNLanczosMethod.cpp:1908
CKNLanczosMethod::EIGENVALUE_RESULT::pEigenVectors
double ** pEigenVectors
Definition: KNLanczosMethod.h:33
CKNMatrixOperation::CKNVector::SetAt
void SetAt(unsigned int nIndex, CKNComplex value)
Set element value in specific index, Call by value.
Definition: KNMatrixOperation.cpp:54
CKNLanczosMethod::EIGENVALUE_RESULT::nEigenValueCount
unsigned int nEigenValueCount
Definition: KNLanczosMethod.h:36
FREE_MEM
#define FREE_MEM(pointer)
Macro for memory allocation and assign null value.
Definition: CKNGlobal.h:20
CKNMatrixOperation::CKNVector::GetSize
unsigned int GetSize()
Return Vector elements size.
Definition: KNMatrixOperation.h:71
CKNLanczosMethod::~CKNLanczosMethod
~CKNLanczosMethod()
Destructor.
Definition: KNLanczosMethod.cpp:67
CKNLanczosMethod::EIGENVALUE_RESULT::nEigenValueCountForMemeory
unsigned int nEigenValueCountForMemeory
Definition: KNLanczosMethod.h:41
CKNLanczosMethod::FinalizeTemporaryArrayAndVector
void FinalizeTemporaryArrayAndVector()
Finalize temporary eigenvalue arrays and vectors.
Definition: KNLanczosMethod.cpp:95
CKNLanczosMethod::ExtractDoubleValues
void ExtractDoubleValues(double *pTarget, double *pSource, unsigned int nSrcCount, int *pFilter, unsigned int nFilterCount, bool bExclusive)
Extract value by condition that described in filter.
Definition: KNLanczosMethod.cpp:1421
ERROR_MALLOC
const unsigned long ERROR_MALLOC
Error code that means error occur during memory allocation.
Definition: CKNGlobal.h:62
CKNComplex
This class for complex operation and saving value.
Definition: KNComplex.h:18
KNMPIManager.h
MPI Mangement class.
CKNMPIManager::BarrierAllComm
static void BarrierAllComm()
Is Multilevel MPI Setting.
Definition: KNMPIManager.cpp:791
CKNLanczosMethod::ReleaseResult
static void ReleaseResult(LPEIGENVALUE_RESULT lpResult, bool bReleaseStruct)
Release memory for lanczos method result.
Definition: KNLanczosMethod.cpp:1480
CKNTimeMeasurement::MeasurementStart
static void MeasurementStart(MEASUREMENT_INDEX index)
Measurement start for part.
Definition: KNTimeMeasurement.cpp:73
CKNLanczosMethod::DistinguishClusterOfEigenvalueEx
int DistinguishClusterOfEigenvalueEx(int nEigenValueCount, double *pEigenValues, double *pEigenVectors, bool *pbValidEigenValues, int nIterationCount)
Distinguish clusters values.
Definition: KNLanczosMethod.cpp:1233
CKNLanczosMethod::LanczosIteration
LPEIGENVALUE_RESULT LanczosIteration()
Doing lanczos basic iteration.
Definition: KNLanczosMethod.cpp:173
CKNMPIManager::IsRootRank
static bool IsRootRank()
Get Total node count.
Definition: KNMPIManager.cpp:217
CKNLanczosMethod::LanczosIterationLoop
void LanczosIterationLoop(LPEIGENVALUE_RESULT lpResult, CKNMatrixOperation::CKNVector *V1, unsigned int nIterationCount, CKNComplex *pAlpha, double *pAlphaReal, double *pBeta, double *pWj, double *pWjm1, double *pWjp1, bool bMakeEigvVector=false)
Doing lanczos basic iteration.
Definition: KNLanczosMethod.cpp:303
CKNLanczosMethod::m_bStop
static bool m_bStop
Determind stop iteration before end of iteration count.
Definition: KNLanczosMethod.h:131
COMMAND_SIZE
#define COMMAND_SIZE
Definition: CKNGlobal.h:102
CKNLanczosMethod::RangeChecking
int RangeChecking(int nEigenValueCount, double *pEigenValues, double *pEiegnVectors, double *pRangeCheckingEigenValues, double *pRangeCheckingVectors, int nIterationCount)
Checking eigenvalue range.
Definition: KNLanczosMethod.cpp:1173
CKNMatrixOperation::IsSame
static bool IsSame(double operand1, double operand2, double tol)
Compare two double variable.
Definition: KNMatrixOperation.cpp:2789
CKNMatrixOperation::Gram_schmidt
static int Gram_schmidt(CKNVector *pVect1, CKNVector *pVect2)
Doing gam schmidt orthogonalization.
Definition: KNMatrixOperation.cpp:2814
CKNLanczosMethod::EigenValueSolver
int EigenValueSolver(unsigned int nIterationCount, double *pAlpha, double *pBeta, double *pEigenValues, double *pEigenVectors)
EigenValue Solving.
Definition: KNLanczosMethod.cpp:1321
CKNMPIManager::IsMultiLevelMPI
static bool IsMultiLevelMPI()
Get MPI_Comm.
Definition: KNMPIManager.h:68
CKNLanczosMethod::RESULT_SORT_DATA::nEigenValueFoundIteration
int nEigenValueFoundIteration
Definition: KNLanczosMethod.h:51
CKNLanczosMethod::ConvergenceChecking
int ConvergenceChecking(int nEigenValueCount, double *pEigenValues, double *pEiegnVectors, double *pConvergedEigenValues, double *pConvergedEigenVectors, double fANorm, double *pBeta, int nIterationCount)
Checking convergence criteria.
Definition: KNLanczosMethod.cpp:1138
DO_NOT_CONVERGENCE_CHECKING
#define DO_NOT_CONVERGENCE_CHECKING
Convergernece checking option default value.
Definition: CKNGlobal.h:46
CKNMatrixOperation::CKNVector
This class for describing vector for Lanczos method.
Definition: KNMatrixOperation.h:60
CKNLanczosMethod::ShowLanczosResult
static void ShowLanczosResult(CKNLanczosMethod::LPEIGENVALUE_RESULT lpResult, bool bCalculateEigenVectors, bool bCalculateWaveFunction, double *pKValue, int nRepeatCount)
Save calculating result into file.
Definition: KNLanczosMethod.cpp:1822
CKNMatrixOperation::VVDot
static bool VVDot(CKNVector *pVector1, CKNVector *pVector2, CKNComplex *pResult)
Between vectors dot product operation.
Definition: KNMatrixOperation.cpp:2691
CKNMPIManager::IsDeflationRoot
static bool IsDeflationRoot()
Checking is root rank of Lanczos computation.
Definition: KNMPIManager.h:72
SEND_BACK_EV
#define SEND_BACK_EV
Definition: CKNGlobal.h:97