source: pcp/src/gramsch.c@ 201832

Last change on this file since 201832 was ef282f, checked in by Frederik Heber <heber@…>, 18 years ago

verbosity changed()
  • Property mode set to 100644
File size: 47.7 KB
Line 
1/** \file gramsch.c
2 * Gram-Schmidt-Orthonormalisation.
3 * Herein are all the functions necessary to orthogonalize and normalize the wave
4 * functions OnePsiElement, such as initialization FirstInitGramSchData(), norm
5 * GramSchNormalize(), scalar product GramSchSP() and the actual Gram-Schmidt-routine
6 * GramSch(). All depending on the current status of the wave function.
7 *
8 Project: ParallelCarParrinello
9 \author Jan Hamaekers
10 \date 2000
11
12 File: gramsch.c
13 $Id: gramsch.c,v 1.70.2.1 2007-04-21 12:49:50 foo Exp $
14*/
15
16#include <stdlib.h>
17#include <stdio.h>
18#include <math.h>
19#include <string.h>
20
21// use double precision fft when we have it
22#ifdef HAVE_CONFIG_H
23#include <config.h>
24#endif
25
26#ifdef HAVE_DFFTW_H
27#include "dfftw.h"
28#else
29#include "fftw.h"
30#endif
31
32#include "data.h"
33#include "errors.h"
34#include "gramsch.h"
35#include "helpers.h"
36#include "myfft.h"
37#include "mymath.h"
38#include "mergesort2.h"
39#include "perturbed.h"
40#include "run.h"
41
42/** Deallocates the defined OnePsiElement datatype.
43 */
44void FreeMPI_OnePsiElement()
45{
46 MPI_Type_free(&MPI_OnePsiElement);
47}
48
49/** Initialization of Gram-Schmidt-Orthogonalization.
50 * \param *P Problem at hand
51 * \param *Psi wave functions
52 * \sa RemoveEverything()
53 */
54void FirstInitGramSchData(struct Problem *P, struct Psis *Psi) {
55 int i, type;
56 int GramSchLocalNo = Psi->LocalNo+1;
57 MPI_Datatype type1[10] = { MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_UB}; // type of each OPE array element
58 int blocklen1[10] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; // block length of each element within the OPE array
59 MPI_Aint base, disp1[10]; // holds adresses in memory
60 struct OnePsiElement OPE[2];
61
62 /// Create MPI_OnePsiElement, simulacrum of OnePsiElement, enabling exchange of these among the processes
63 // store adresses of its various elements in disp1 array
64 MPI_Address( OPE, disp1);
65 MPI_Address( &OPE[0].my_color_comm_ST_Psi, disp1+1);
66 MPI_Address( &OPE[0].MyLocalNo, disp1+2);
67 MPI_Address( &OPE[0].MyGlobalNo, disp1+3);
68 MPI_Address( &OPE[0].PsiGramSchStatus, disp1+4);
69 MPI_Address( &OPE[0].PsiType, disp1+5);
70 MPI_Address( &OPE[0].PsiFactor, disp1+6);
71 MPI_Address( &OPE[0].PsiReciNorm2, disp1+7);
72 MPI_Address( &OPE[0].DoBrent, disp1+8);
73 MPI_Address( OPE+1, disp1+9);
74 base = disp1[0];
75 for (i=0; i < 10; i++) disp1[i] -= base; // make the adresses of OPE elements relativ to base -> byte displacement of each entry
76 MPI_Type_struct( 10, blocklen1, disp1, type1, &MPI_OnePsiElement); // creates MPI_OnePsiElement as an MPI_struct(ure)
77 MPI_Type_commit( &MPI_OnePsiElement); // commits new data type, now it's usable
78 if (P->Call.out[NormalOut]) fprintf(stderr, "(%i)FirstInitGramSchData\n", P->Par.me);
79
80 /// Allocates and fills Psis::AllLocalNo (MPI_Allgathered from all other processes).
81 Psi->AllLocalNo = (int *)
82 Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"FirstInitGramSchData: Psi->AllLocalNo");
83 MPI_Allgather ( &GramSchLocalNo, 1, MPI_INT,
84 Psi->AllLocalNo, 1, MPI_INT, P->Par.comm_ST_PsiT );
85
86 /// Calculates from this Psis::MaxPsiOfType.
87 Psi->MaxPsiOfType = 0;
88 for (i=0;i<P->Par.Max_me_comm_ST_PsiT;i++) Psi->MaxPsiOfType += Psi->AllLocalNo[i]-1; // sum up all local (orthogonalizable) Psis in the transposed communicator PsiT
89 if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) MaxPsiOfType = %i\n",P->Par.me, Psi->MaxPsiOfType);
90
91 /// Calculates from this Psis::MaxPsiOfType and at which index this process' Psis start Psis::MyStartNo.
92 Psi->MyStartNo = 0;
93 for (i=0;i<P->Par.me_comm_ST_PsiT;i++) Psi->MyStartNo += Psi->AllLocalNo[i]; // where do my Psis start
94 if (P->Call.out[NormalOut]) fprintf(stderr,"(%i) MyStartNo = %i\n",P->Par.me, Psi->MyStartNo);
95
96 //fprintf(stderr,"(%i) OtherPsiLocalNo %d\n",P->Par.me, RecvCount);
97 /// Allocates arrays Psis::AllPsiStatus, Psis::AllPsiStatusForSort and Psis::LocalPsiStatus (up 'til Extra in PsiTagType)
98 Psi->AllPsiStatus = (struct OnePsiElement *)
99 Malloc(sizeof(struct OnePsiElement)*(Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT),"FirstInitGramSchData: Psi->AllPsiStatus");
100 Psi->AllPsiStatusForSort = (struct OnePsiElement *)
101 Malloc(sizeof(struct OnePsiElement)*(Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT+1),"FirstInitGramSchData: Psi->AllPsiStatusForSort");
102 Psi->LocalPsiStatus = (struct OnePsiElement *)
103 Malloc(sizeof(struct OnePsiElement)*GramSchLocalNo,"FirstInitGramSchData: Psi->LocalPsiStatus");
104 /// Psis::LocalPsiStatus is initialized and distributed among all processes as Psis::AllPsiStatus.
105 for (i=0;i<GramSchLocalNo;i++) {
106 Psi->LocalPsiStatus[i].me_comm_ST_Psi = P->Par.me_comm_ST_Psi;
107 Psi->LocalPsiStatus[i].my_color_comm_ST_Psi = P->Par.my_color_comm_ST_Psi;
108 Psi->LocalPsiStatus[i].MyLocalNo = i;
109 Psi->LocalPsiStatus[i].MyGlobalNo = Psi->MyStartNo + i;
110 Psi->LocalPsiStatus[i].DoBrent = 4;
111 switch (Psi->PsiST) { // set occupation number for the regular local, one extra(!) per process (without current one!) and the additional orbitals (the "latterest" ;) are set to zero of course) (NOTE: extra orbit must always be the very last one (that's why Par->.. - 1)
112 case SpinDouble:
113 for (type=Occupied;type<=Extra;type++) {
114 if (i >= Psi->TypeStartIndex[type] && i < Psi->TypeStartIndex[type+1]) {
115 Psi->LocalPsiStatus[i].PsiType = type;
116 Psi->LocalPsiStatus[i].PsiGramSchStatus = (int)(type != Occupied ? NotUsedToOrtho : NotOrthogonal); // extra or occupied wave function
117 }
118 }
119 if (Psi->LocalPsiStatus[i].PsiType != UnOccupied)
120 Psi->LocalPsiStatus[i].PsiFactor = 2.0;
121 else
122 Psi->LocalPsiStatus[i].PsiFactor = 1.0;
123 break;
124 case SpinUp:
125 for (type=Occupied;type<=Extra;type++) {
126 if (i >= Psi->TypeStartIndex[type] && i < Psi->TypeStartIndex[type+1]) {
127 Psi->LocalPsiStatus[i].PsiType = type;
128 Psi->LocalPsiStatus[i].PsiGramSchStatus = (int)(type != Occupied ? NotUsedToOrtho : NotOrthogonal); // extra or occupied wave function
129 }
130 }
131 Psi->LocalPsiStatus[i].PsiFactor = 1.0;
132 break;
133 case SpinDown:
134 for (type=Occupied;type<=Extra;type++) {
135 if (i >= Psi->TypeStartIndex[type] && i < Psi->TypeStartIndex[type+1]) {
136 Psi->LocalPsiStatus[i].PsiType = type;
137 Psi->LocalPsiStatus[i].PsiGramSchStatus = (int)(type != Occupied ? NotUsedToOrtho : NotOrthogonal); // extra or occupied wave function
138 }
139 }
140 Psi->LocalPsiStatus[i].PsiFactor = 1.0;
141 break;
142 }
143 Psi->LocalPsiStatus[i].PsiReciNorm2 = 0.0;
144 }
145
146 // Update AllPsiStatus from changed LocalPsiStatus
147 UpdateGramSchAllPsiStatus(P,Psi);
148
149 /// Psis::TempSendA, Psis::AllActualLocalPsiNo and Psis::AllOldActualLocalPsiNo are allocated, the latter two zeroed.
150 Psi->TempSendA = (int *)
151 Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"FirstInitGramSchData: Psi->TempSendA");
152 Psi->AllActualLocalPsiNo = (int *)
153 Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"FirstInitGramSchData: Psi->AllActualLocalPsiNo");
154 Psi->AllOldActualLocalPsiNo = (int *)
155 Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"FirstInitGramSchData: Psi->AllOldActualLocalPsiNo");
156 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
157 Psi->AllActualLocalPsiNo[i] = 0;
158 Psi->AllOldActualLocalPsiNo[i] = 0;
159 }
160}
161
162/** Normalize the coefficients of a given wave function.
163 * Calculates the norm (see GramSchGetNorm2()) and divides each (for all reciprocal grid
164 * vectors) complex coefficient by the norm.
165 * \param *P Problem at hand
166 * \param *Lev LatticeLevel structure
167 * \param *LPsiDat Array of complex wave function coefficients
168 * \param PsiSP If norm already calculated, can be passed on here, otherweise (== 0.0) is calculated
169 * \return Squared norm of wave function
170 */
171double GramSchNormalize(const struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDat, double PsiSP) {
172 double LocalSP=0.0;
173 int i,s = 0;
174 /* Falls PsiSP == 0.0 dann noch SP berechnen */
175 if (PsiSP == 0.0) { // see GramSchGetNorm2()
176 if (Lev->GArray[0].GSq == 0.0) {
177 LocalSP += LPsiDat[0].re*LPsiDat[0].re;
178 s++;
179 }
180 for (i=s; i < Lev->MaxG; i++) {
181 LocalSP += 2*(LPsiDat[i].re*LPsiDat[i].re+LPsiDat[i].im*LPsiDat[i].im);
182 }
183 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
184 }
185 if ((PsiSP < MYEPSILON) && (P->Call.out[PsiOut])) fprintf(stderr,"GramSchNormalize: PsiSP = %lg\n",PsiSP);
186 PsiSP = sqrt(PsiSP); // take square root
187 for (i=0; i < Lev->MaxG; i++) { // and divide each coefficient by the norm
188 LPsiDat[i].re /= PsiSP;
189 LPsiDat[i].im /= PsiSP;
190 }
191 return(PsiSP);
192}
193
194/** Calculate squared norm of given wave function coefficients.
195 * Go through each node of the reciprocal vector grid, calculate the complex product for this
196 * coefficient and sum up, gathering the results from all processes before return - remember
197 * that the coefficients are - for the parallel calculation of the fft - split up among the
198 * processes.
199 * \param *P Problem at hand
200 * \param *Lev LatticeLevel structure
201 * \param *LPsiDat array over G of complex i-th wave function coefficients \f$c_{i,G}\f$
202 * \return \f$\sum_G c_{i,G} /cdot {c_{i,G}}^{\ast}\f$
203 */
204double GramSchGetNorm2(const struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDat) {
205 double LocalSP=0.0, PsiSP=0.0;
206 int i,s = 0;
207 /* Falls PsiSP == 0.0 dann noch SP berechnen */
208 if (LPsiDat != NULL) {
209 if (Lev->GArray[0].GSq == 0.0) {
210 LocalSP += LPsiDat[0].re*LPsiDat[0].re;
211 s++;
212 }
213 for (i=s; i < Lev->MaxG; i++) {
214 LocalSP += 2*(LPsiDat[i].re*LPsiDat[i].re+LPsiDat[i].im*LPsiDat[i].im);
215 }
216 // send local result to all processes and received summed from all into PsiSP
217 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
218 }
219 return(PsiSP);
220}
221
222/** Scalar Product of two arrays of wave function coefficients.
223 * Goes through each reciprocal grid vectors and calculates the complex product
224 * between the two coefficients, summing up, MPI_Allreducing and returning.
225 * (See also GramSchGetNorm2())
226 * \param *P Problem at hand
227 * \param *Lev LatticeLevel structure
228 * \param *LPsiDatA first array of wave function coefficients
229 * \param *LPsiDatB second array of wave function coefficients
230 * \return \f$\sum_G c_{a,G} \cdot c_{b,G}^{\ast}\f$
231 */
232static double GramSchSP(const struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDatA, fftw_complex *LPsiDatB) {
233 double LocalSP=0.0,PsiSP;
234 int i,s = 0;
235 if (Lev->GArray[0].GSq == 0.0) {
236 LocalSP += LPsiDatA[0].re*LPsiDatB[0].re;
237 s++;
238 }
239 for (i=s; i < Lev->MaxG; i++) { // go through all nodes and calculate complex scalar product
240 LocalSP += 2*(LPsiDatA[i].re*LPsiDatB[i].re+LPsiDatA[i].im*LPsiDatB[i].im);
241 }
242 // send local result to all processes and received summed from all into PsiSP
243 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
244 return(PsiSP);
245}
246
247/** Sort criteria for natueralmergesort(): Returns re-ordered OnePsiElement::PsiGramSchStatus.
248 * The current status in the Gram-Schmidt-Orthonormalization is returned as sort criteria.
249 * \param *a OnePsiElement
250 * \param i i-th wave function
251 * \param *Args
252 * \return integer value for each PsiGramSchStatusType, from IsOrthonormal (0) up to NotOrthogonal(2) and NotUsedToOrtho(3)
253 * \note The enum PsiGramSchStatusType is not simply copied due to a different ordering in the enumeration other than used here.
254 */
255static double GetKeyOnePsi(void *a, int i, void *Args) {
256 double res=-1;
257 switch ((enum PsiGramSchStatusType)((struct OnePsiElement *)a)[i].PsiGramSchStatus) {
258 case NotOrthogonal:
259 res = 2.;
260 break;
261 case IsOrthogonal:
262 res = 1.;
263 break;
264 case IsOrthonormal:
265 res = 0.;
266 break;
267 case NotUsedToOrtho:
268 res = 100.; // extra before unoccupied and perturbed ones
269 break;
270 }
271 switch (((struct OnePsiElement *)a)[i].PsiType) {
272 case Occupied:
273 res += 0.;
274 break;
275 case UnOccupied:
276 res += 10.;
277 break;
278 case Perturbed_P0:
279 case Perturbed_P1:
280 case Perturbed_P2:
281 case Perturbed_RxP0:
282 case Perturbed_RxP1:
283 case Perturbed_RxP2:
284 res += 20.;
285 break;
286 case Extra:
287 res += 30.;
288 break;
289 }
290 return(res);
291}
292
293/** Sort criteria for natueralmergesort(): Returns the global number of the Psi among all.
294 * \param *a OnePsiElement
295 * \param i i-th wave function
296 * \param *Args unused, for contingency with GetKeyOnePsi()
297 * \return \a i-th OnePsiElement::MyGlobalNo
298 */
299static double GetKeyOnePsi2(void *a, int i, void *Args) {
300 return(((struct OnePsiElement *)a)[i].MyGlobalNo);
301}
302
303/** Copies wave function OnePsiElement.
304 * Copy each entry in OnePsiElement structure from \a b[j] to \a a[i].
305 * \param *a destination OnePsiElement array
306 * \param i i-th element to be overwritten
307 * \param *b source OnePsiElement array
308 * \param j j-th element's entries to be copied
309 */
310static void CopyElementOnePsi(void *a, int i, void *b, int j)
311{
312 ((struct OnePsiElement *)a)[i].me_comm_ST_Psi = ((struct OnePsiElement *)b)[j].me_comm_ST_Psi;
313 ((struct OnePsiElement *)a)[i].my_color_comm_ST_Psi = ((struct OnePsiElement *)b)[j].my_color_comm_ST_Psi;
314 ((struct OnePsiElement *)a)[i].MyLocalNo = ((struct OnePsiElement *)b)[j].MyLocalNo;
315 ((struct OnePsiElement *)a)[i].MyGlobalNo = ((struct OnePsiElement *)b)[j].MyGlobalNo;
316 ((struct OnePsiElement *)a)[i].PsiGramSchStatus = ((struct OnePsiElement *)b)[j].PsiGramSchStatus;
317 ((struct OnePsiElement *)a)[i].PsiType = ((struct OnePsiElement *)b)[j].PsiType;
318 ((struct OnePsiElement *)a)[i].PsiFactor = ((struct OnePsiElement *)b)[j].PsiFactor;
319 ((struct OnePsiElement *)a)[i].PsiReciNorm2 = ((struct OnePsiElement *)b)[j].PsiReciNorm2;
320}
321
322
323/** Performs Gram-Schmidt-Orthonormalization on all Psis.
324 * Herein the known Gram-Schmidt-Orthogonalization (with subsequent normalization) is implemented in a
325 * parallel way. The problem arises due to the fact that the complex wave function coefficients are not
326 * all accessible from one process, but are shared among them. Thus there are four different cases to
327 * deal with - where O is one orthogonal Psi and P the Psi currently to be orthogonalized:
328 * -# O and P are local\n
329 * The projection is simply calculated via scalar product and subtracted from P.
330 * -# O is local, P not\n
331 * P is received from the respective process and the projetion calculated, noting down this
332 * value for later sending it back to this respective process owning the P coefficients,
333 * who will substract them
334 * -# O is not local, however P is\n
335 * Send the coefficient to every process in need of them and in the end gather projections to
336 * be subtracted from our local P.
337 * -# O and P are not local\n
338 * Nothing to do.
339 *
340 * Afterwards, a division by the norm of the Psi may additionally be called in for. The current status of
341 * a Psi is always noted in OnePsiElement::PsiGramSchStatus.
342 * \param *P Problem at hand
343 * \param *Lev LatticeLevel structure
344 * \param *Psi wave functions structure Psis
345 * \param ToDo states what to do in this function: Orthogonalize or Orthonormalize
346 */
347void GramSch(struct Problem *P, struct LatticeLevel *Lev, struct Psis *Psi, enum PsiGramSchToDoType ToDo)
348{
349 int i, j, k, TempRecv, TempSend, RecvSource;
350 //int ResetNo=0;
351 double GlobalSP;
352 struct RunStruct *R = &P->R;
353 struct OnePsiElement *OnePsi = NULL, *LOnePsi = NULL, *ROnePsi = NULL, *RLOnePsi = NULL;
354 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
355 fftw_complex *Temp, *Temp2;
356 int *TempSendA = Psi->TempSendA;
357 MPI_Status status;
358 // Mergesort the wavefunction by their current status from 0 to all plus all extra ones (one for each process)
359 naturalmergesort(Psi->AllPsiStatus,Psi->AllPsiStatusForSort,0,Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT-1,&GetKeyOnePsi,NULL,&CopyElementOnePsi);
360 //fprintf(stderr,"(%i) GramSch: ",P->Par.me);
361 for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // then go through each of the ToDo-order sorted Psis (Each Psi plus an extra one from each process)
362 OnePsi = &Psi->AllPsiStatus[i]; // Mark OnePsi wave function
363 if (OnePsi->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // stored in this process? => L means local
364 LOnePsi = &Psi->LocalPsiStatus[OnePsi->MyLocalNo];
365 else
366 LOnePsi = NULL;
367
368 switch ((enum PsiGramSchStatusType)OnePsi->PsiGramSchStatus) { // depending on their ToDo-status do ...
369 case NotOrthogonal: // ORTHOGONALIZE!
370 //fprintf(stderr,"(%i) ", Psi->AllPsiStatus[i].MyGlobalNo);
371 //fprintf(stderr,"Orthogonalizing %i, status was: (L)%i\t(A)%i!\n", OnePsi->MyGlobalNo, Psi->LocalPsiStatus[OnePsi->MyLocalNo].PsiGramSchStatus, OnePsi->PsiGramSchStatus);
372 if (LOnePsi != NULL) { // if current Psi is local, copy (reciprocal) complex coefficients from LocalPsi to TempPsi
373 memcpy(Lev->LPsi->TempPsi, Lev->LPsi->LocalPsi[OnePsi->MyLocalNo], ElementSize*Lev->MaxG*sizeof(double));
374 }
375 Temp = Lev->LPsi->TempPsi2; // another complex coefficients array (reciprocal) ...
376 SetArrayToDouble0((double *)Temp, Lev->MaxG*2); // ... which is zeroed
377 TempRecv = 0; // count how often a needed local current Psi has been received (and thus if it has been already)
378 TempSend = 0; // count how often a local current Psi has been sent to other processes
379 for(k=0; k < P->Par.Max_me_comm_ST_PsiT; k++) TempSendA[k] = 0; // zero array counting how often a process has sent its local Psi to others
380
381 for (j=i-1; j >= 0; j--) { // go through all wave functions from the one before the current downto first one
382 ROnePsi = &Psi->AllPsiStatus[j]; // get the Psi that should be orthogonal to it
383 if (ROnePsi->PsiType <= UnOccupied) { // only orthogonalize against non-perturbed wave functions
384 if (ROnePsi->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // stored in this process? => L means local
385 RLOnePsi = &Psi->LocalPsiStatus[ROnePsi->MyLocalNo];
386 else
387 RLOnePsi = NULL;
388// if (((OnePsi->PsiType == Extra && (R->CurrentMin <= UnOccupied || ((LOnePsi != NULL && RLOnePsi != NULL) && ROnePsi->MyLocalNo == R->ActualLocalPsiNo)))
389// || OnePsi->PsiType <= UnOccupied)) {
390 if ((ROnePsi->PsiType == Occupied)
391 || ((ROnePsi->PsiType == UnOccupied) && (OnePsi->PsiType == UnOccupied || OnePsi->PsiType == Extra))
392 || ((LOnePsi != NULL && RLOnePsi != NULL) && ROnePsi->MyLocalNo == R->ActualLocalPsiNo)) {
393 // occupied are orthogonal to occupied
394 // unoccupied are orthogonal to occupied and unoccupied
395 // perturbed are orthogonal to occupied, unoccupied and to their (process-) specific extra
396 // extra is orthogonal dependent on R->CurrentMin (to occupied, occupied&unoccupied, occupied&specific perturbed)
397 if (RLOnePsi != NULL && LOnePsi != NULL) { // if both are stored locally in this process
398 GlobalSP = GramSchSP(P,Lev,Lev->LPsi->LocalPsi[ROnePsi->MyLocalNo],Lev->LPsi->LocalPsi[OnePsi->MyLocalNo]); // scalar product of the two
399 GlobalSP *= RLOnePsi->PsiReciNorm2; // divide by norm
400 Temp = Lev->LPsi->TempPsi;
401 Temp2 = Lev->LPsi->LocalPsi[ROnePsi->MyLocalNo];
402 for(k=0; k < Lev->MaxG; k++) { // orthogonalize it (subtract the projected part, real and imaginary)
403 Temp[k].re -= GlobalSP*Temp2[k].re;
404 Temp[k].im -= GlobalSP*Temp2[k].im;
405 } // the orthogonalized wave function of LocalPsi resides now in Temp!
406 }
407 if (RLOnePsi != NULL && LOnePsi == NULL) { // if the current Psi is not local, the one to which it ought be orthogonal however is
408 /* Recv i and put it to jLocal in TempPsi */
409 if (TempRecv == 0) {
410 MPI_Recv( Lev->LPsi->TempPsi, Lev->MaxG*ElementSize, MPI_DOUBLE, OnePsi->my_color_comm_ST_Psi, GramSchTag1, P->Par.comm_ST_PsiT, &status );
411 TempRecv++;
412 }
413 GlobalSP = GramSchSP(P,Lev,Lev->LPsi->LocalPsi[ROnePsi->MyLocalNo],Lev->LPsi->TempPsi);
414 GlobalSP *= RLOnePsi->PsiReciNorm2;
415 Temp = Lev->LPsi->TempPsi2;
416 Temp2 = Lev->LPsi->LocalPsi[ROnePsi->MyLocalNo];
417 for(k=0; k < Lev->MaxG; k++) {
418 Temp[k].re -= GlobalSP*Temp2[k].re;
419 Temp[k].im -= GlobalSP*Temp2[k].im;
420 } // the negative orthogonal projection resides in Temp (not the local wave function part!)
421 }
422 if (RLOnePsi == NULL && LOnePsi != NULL) { // if the current Psi is local, the one to which it ought be orthogonal yet is not
423 /* Send i to jLocal in TempPsi */
424 if (TempSendA[ROnePsi->my_color_comm_ST_Psi] == 0) { // just send it out to everyone who needs it
425 MPI_Send( Lev->LPsi->LocalPsi[OnePsi->MyLocalNo], Lev->MaxG*ElementSize, MPI_DOUBLE, ROnePsi->my_color_comm_ST_Psi, GramSchTag1, P->Par.comm_ST_PsiT);
426 TempSend++;
427 TempSendA[ROnePsi->my_color_comm_ST_Psi]++; // note that coefficients were sent once more to this process
428 }
429 }
430 }
431 }
432 }
433 if (LOnePsi != NULL) { // holds the current local Psi (TempPsi) and receives results from all other which "ought be orthogonal to this one"
434 /* Hat was in TempPsi und ist local*/
435 for (j=0; j < TempSend; j++) { // each of the recipients before should send something back now
436 MPI_Probe( MPI_ANY_SOURCE, GramSchTag2, P->Par.comm_ST_PsiT, &status );
437 RecvSource = status.MPI_SOURCE;
438 MPI_Recv( Lev->LPsi->TempPsi2, Lev->MaxG*ElementSize, MPI_DOUBLE, RecvSource, GramSchTag2, P->Par.comm_ST_PsiT, &status );
439 Temp2 = Lev->LPsi->TempPsi2;
440 Temp = Lev->LPsi->TempPsi;
441 for(k=0; k < Lev->MaxG; k++) { // sum received projetion onto (temporary) local wave function
442 Temp[k].re += Temp2[k].re;
443 Temp[k].im += Temp2[k].im;
444 }
445 }
446 Temp2 = Lev->LPsi->TempPsi;
447 Temp = Lev->LPsi->LocalPsi[OnePsi->MyLocalNo];
448 memcpy(Temp, Temp2, sizeof(fftw_complex)*Lev->MaxG); // finally copy back onto original one
449 }
450 if (LOnePsi == NULL && TempRecv) { // has calculated a projection to another Psi (TempPsi2) and sends it to the respective (local) owner of the current one
451 /* Hat was in TempPsi2 und ist nicht local*/
452 MPI_Send( Lev->LPsi->TempPsi2, Lev->MaxG*ElementSize, MPI_DOUBLE, OnePsi->my_color_comm_ST_Psi, GramSchTag2, P->Par.comm_ST_PsiT);
453 }
454 /*if (LOnePsi != NULL) { // finally we set the status of our local (multi-projection subtracted) Psi
455 //fprintf(stderr,"Setting L-Status of %i to %i\n",LOnePsi->MyGlobalNo, IsOrthogonal);
456 LOnePsi->PsiGramSchStatus = (int)IsOrthogonal;
457 }
458 fprintf(stderr,"Setting A-Status of %i to %i\n",OnePsi->MyGlobalNo, IsOrthogonal);
459 OnePsi->PsiGramSchStatus = (int)IsOrthogonal; // and also set the status in all processes for this Psi*/
460 // note: There is no break here, normalization will be performed right away!
461 //fprintf(stderr,"-> ");
462 case IsOrthogonal: // NORMALIZE!
463 //fprintf(stderr,"%i ", Psi->AllPsiStatus[i].MyGlobalNo);
464 switch (ToDo) {
465 case Orthonormalize: // ... normalize and store 1 as norm
466 if (LOnePsi != NULL) {
467 //fprintf(stderr,"Setting L-Status of %i to %i\n",LOnePsi->MyLocalNo, IsOrthonormal);
468 LOnePsi->PsiGramSchStatus = (int)IsOrthonormal; // set status and ...
469 GramSchNormalize(P,Lev,Lev->LPsi->LocalPsi[OnePsi->MyLocalNo],0.0); // ... do it
470 /*LOnePsi->PsiReciNorm2 = GramSchGetNorm2(P,Lev,Lev->LPsi->LocalPsi[OnePsi->MyLocalNo]);
471 LOnePsi->PsiReciNorm2 = 1./LOnePsi->PsiReciNorm2;*/
472 LOnePsi->PsiReciNorm2 = 1.;
473 }
474 //fprintf(stderr,"Setting A-Status of %i to %i\n",OnePsi->MyGlobalNo, IsOrthonormal);
475 OnePsi->PsiGramSchStatus = (int)IsOrthonormal;
476 break;
477 case Orthogonalize: // ... calculate norm and store
478 if (LOnePsi != NULL) {
479 //fprintf(stderr,"Setting L-Status of %i to %i\n",LOnePsi->MyLocalNo, IsOrthogonal);
480 LOnePsi->PsiGramSchStatus = (int)IsOrthogonal;
481 LOnePsi->PsiReciNorm2 = GramSchGetNorm2(P,Lev,Lev->LPsi->LocalPsi[OnePsi->MyLocalNo]);
482 if ((LOnePsi->PsiReciNorm2 < MYEPSILON) && (P->Call.out[PsiOut])) fprintf(stderr,"GramSch: LOnePsi->PsiReciNorm2 Nr. %d = %lg\n",LOnePsi->MyGlobalNo,LOnePsi->PsiReciNorm2);
483 LOnePsi->PsiReciNorm2 = 1./LOnePsi->PsiReciNorm2;
484 }
485 //fprintf(stderr,"Setting A-Status of %i to %i\n",OnePsi->MyGlobalNo, IsOrthogonal);
486 OnePsi->PsiGramSchStatus = (int)IsOrthogonal;
487 break;
488 }
489 break;
490 case IsOrthonormal: // NOTHING TO DO ANY MORE!
491 //fprintf(stderr,"%i ", Psi->AllPsiStatus[i].MyGlobalNo);
492 break;
493 case NotUsedToOrtho:
494 //fprintf(stderr,"[%i] ", Psi->AllPsiStatus[i].MyGlobalNo);
495 break;
496 }
497 }
498 //fprintf(stderr,"\n");
499 /* Reset */
500 naturalmergesort(Psi->AllPsiStatus,Psi->AllPsiStatusForSort,0,Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT-1,&GetKeyOnePsi2,NULL,&CopyElementOnePsi);
501/* fprintf(stderr,"Setting L-Status of %i to %i\n",Psi->LocalNo, NotUsedToOrtho);
502 Psi->LocalPsiStatus[Psi->LocalNo].PsiGramSchStatus = (int)NotUsedToOrtho;
503 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
504 ResetNo += Psi->AllLocalNo[i];
505 OnePsi = &Psi->AllPsiStatus[ResetNo-1];
506 fprintf(stderr,"Setting A-Status of %i to %i\n",OnePsi->MyGlobalNo, NotUsedToOrtho);
507 OnePsi->PsiGramSchStatus = (int)NotUsedToOrtho;
508 }*/
509}
510
511/** Reset status of Gram-Schmidt-Orthogonalization for each and every Psi.
512 * Sets all locally accessible Psis::LocalPsiStatus to PsiGramSchStatusType::NotOrthogonal
513 * and the norm to zero, except the last (extra) and unoccupied ones which are NotUsedToOrtho, then
514 * do the same for all Psis::AllPsiStatus (again exception for extra and unoccupied ones).
515 * \param *P Problem at hand
516 * \param *Psi wave functions structure Psis
517 */
518void ResetGramSch(const struct Problem *P, struct Psis *Psi)
519{
520 int i,j, ResetNo=0;
521 struct OnePsiElement *OnePsi = NULL;
522 for (i=0; i < Psi->LocalNo; i++) { // go through all local Psis
523 Psi->LocalPsiStatus[i].PsiGramSchStatus = (Psi->LocalPsiStatus[i].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
524 //fprintf(stderr,"Setting L-Status of %i to %i\n",i, Psi->LocalPsiStatus[i].PsiGramSchStatus);
525 Psi->LocalPsiStatus[i].PsiReciNorm2 = 0.0;
526 }
527 //fprintf(stderr,"Setting L-Status of %i to %i\n",Psi->LocalNo, NotUsedToOrtho);
528 Psi->LocalPsiStatus[Psi->LocalNo].PsiGramSchStatus = (int)NotUsedToOrtho; // extra wave function
529 Psi->LocalPsiStatus[Psi->LocalNo].PsiReciNorm2 = 0.0;
530 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
531 for (j=ResetNo; j < ResetNo+Psi->AllLocalNo[i]-1; j++) {
532 Psi->AllPsiStatus[j].PsiGramSchStatus = (Psi->AllPsiStatus[j].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
533 //fprintf(stderr,"Setting A-Status of %i to %i\n",j, Psi->AllPsiStatus[j].PsiGramSchStatus);
534 Psi->AllPsiStatus[j].PsiReciNorm2 = 0.0;
535 }
536 ResetNo += Psi->AllLocalNo[i];
537 OnePsi = &Psi->AllPsiStatus[ResetNo-1];
538 //fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo-1, NotUsedToOrtho);
539 OnePsi->PsiGramSchStatus = (int)NotUsedToOrtho; // extra wave function
540 OnePsi->PsiReciNorm2 = 0.0;
541 }
542}
543
544/** Reset status of Gram-Schmidt-Orthogonalization for each Psi of PsiTagType \a type.
545 * Sets all locally accessible Psis::LocalPsiStatus to PsiGramSchStatusType::NotOrthogonal
546 * and the norm to zero, except the last (extra) and unoccupied ones which are NotUsedToOrtho, then
547 * do the same for all Psis::AllPsiStatus (again exception for extra and unoccupied ones).
548 * \param *P Problem at hand
549 * \param *Psi wave functions structure Psis
550 * \param type PsiTagType of orbitals whose PsiGramSchStatus is to be reset
551 * \param ToDo - set PsiGramSchToDoType for the \a type states to this
552 * \sa ResetGramSch() - same procedure for occupied states
553 */
554void ResetGramSchTagType(const struct Problem *P, struct Psis *Psi, enum PsiTypeTag type, enum PsiGramSchStatusType ToDo)
555{
556 int i,j, ResetNo=0;
557 struct OnePsiElement *OnePsi = NULL;
558 for (i=0; i < Psi->LocalNo; i++) { // go through all local Psis
559 if (Psi->LocalPsiStatus[i].PsiType == type) {
560 //fprintf(stderr,"Setting L-Status of %i to %i\n",i, ToDo);
561 Psi->LocalPsiStatus[i].PsiGramSchStatus = ToDo;
562 switch(ToDo) {
563 case NotOrthogonal:
564 Psi->LocalPsiStatus[i].PsiReciNorm2 = 0.0;
565 break;
566 default:
567 break;
568 }
569 }
570 }
571 //fprintf(stderr,"Setting L-Status of %i to %i\n",Psi->LocalNo, NotUsedToOrtho);
572 Psi->LocalPsiStatus[Psi->LocalNo].PsiGramSchStatus = (int)NotUsedToOrtho; // extra wave function
573 Psi->LocalPsiStatus[Psi->LocalNo].PsiReciNorm2 = 0.0;
574 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
575 for (j=ResetNo; j < ResetNo+Psi->AllLocalNo[i]-1; j++) {
576 if (Psi->AllPsiStatus[j].PsiType == type) {
577 //fprintf(stderr,"Setting A-Status of %i to %i\n",j, ToDo);
578 Psi->AllPsiStatus[j].PsiGramSchStatus = ToDo;
579 switch(ToDo) {
580 case NotOrthogonal:
581 Psi->AllPsiStatus[j].PsiReciNorm2 = 0.0;
582 break;
583 default:
584 break;
585 }
586 }
587 }
588 ResetNo += Psi->AllLocalNo[i];
589 OnePsi = &Psi->AllPsiStatus[ResetNo-1];
590 //fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo-1, NotUsedToOrtho);
591 OnePsi->PsiGramSchStatus = (int)NotUsedToOrtho; // extra wave function
592// OnePsi->PsiReciNorm2 = 0.0;
593 }
594}
595/** Set Gram-Schmidt status of the extra Psis::LocalPsiStatus and Psis::AllPsiStatus Psis to \a PsiGramSchStatus.
596 * The number of these "extra" Psis is Psis::AllLocalNo - 1 for each process.
597 * \param *P Problem at hand
598 * \param *Psi wave functions structure Psis
599 * \param PsiGramSchStatus status to be set
600 */
601void SetGramSchExtraPsi(const struct Problem *P, struct Psis *Psi, enum PsiGramSchStatusType PsiGramSchStatus)
602{
603 int i, ResetNo=0; // offset to the respective local Psis
604 struct OnePsiElement *OnePsi = NULL;
605 //fprintf(stderr,"Setting L-Status of %i to %i\n",Psi->LocalNo, PsiGramSchStatus);
606 Psi->LocalPsiStatus[Psi->LocalNo].PsiGramSchStatus = (int)PsiGramSchStatus;
607 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
608 ResetNo += Psi->AllLocalNo[i];
609 OnePsi = &Psi->AllPsiStatus[ResetNo-1];
610 //fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo-1, PsiGramSchStatus);
611 OnePsi->PsiGramSchStatus = (int)PsiGramSchStatus;
612 }
613}
614
615/** Set Gram-Schmidt status of the actual Psis::LocalPsiStatus and Psis::AllPsiStatus Psis to \a PsiGramSchStatus.
616 * The number of these "extra" Psis is Psis::AllActualLocalPsiNo for each process.
617 * \param *P Problem at hand
618 * \param *Psi wave functions structure Psis
619 * \param PsiGramSchStatus status to be set
620 */
621void SetGramSchActualPsi(const struct Problem *P, struct Psis *Psi, enum PsiGramSchStatusType PsiGramSchStatus)
622{
623 int i, ResetNo=0; // offset to the respective local Psis
624 struct OnePsiElement *OnePsi = NULL;
625 //fprintf(stderr,"Setting L-Status of %i to %i\n",P->R.ActualLocalPsiNo, PsiGramSchStatus);
626 //BUG: Psi->LocalPsiStatus[Psi->LocalNo].PsiGramSchStatus = (int)PsiGramSchStatus;
627 Psi->LocalPsiStatus[P->R.ActualLocalPsiNo].PsiGramSchStatus = (int)PsiGramSchStatus;
628 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
629 OnePsi = &Psi->AllPsiStatus[ResetNo+Psi->AllActualLocalPsiNo[i]];
630 //fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo+Psi->AllActualLocalPsiNo[i], PsiGramSchStatus);
631 OnePsi->PsiGramSchStatus = (int)PsiGramSchStatus;
632 ResetNo += Psi->AllLocalNo[i];
633 }
634}
635
636/** Set Gram-Schmidt status of the former actual Psis::LocalPsiStatus and Psis::AllPsiStatus Psis to \a PsiGramSchStatus.
637 * The former number of these "extra" Psis is Psis::AllOldActualLocalPsiNo for each process.
638 * \param *P Problem at hand
639 * \param *Psi wave functions structure Psis
640 * \param PsiGramSchStatus status to be set
641 */
642void SetGramSchOldActualPsi(const struct Problem *P, struct Psis *Psi, enum PsiGramSchStatusType PsiGramSchStatus)
643{
644 int i, ResetNo=0;
645 struct OnePsiElement *OnePsi = NULL;
646 //fprintf(stderr,"Setting L-Status of %i to %i\n",P->R.OldActualLocalPsiNo, PsiGramSchStatus);
647 Psi->LocalPsiStatus[P->R.OldActualLocalPsiNo].PsiGramSchStatus = (int)PsiGramSchStatus;
648 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
649 OnePsi = &Psi->AllPsiStatus[ResetNo+Psi->AllOldActualLocalPsiNo[i]];
650 //fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo+Psi->AllOldActualLocalPsiNo[i], PsiGramSchStatus);
651 OnePsi->PsiGramSchStatus = (int)PsiGramSchStatus;
652 ResetNo += Psi->AllLocalNo[i];
653 }
654}
655
656/** Updates array Psis::AllActualLocalPsiNo from the RunStruct::ActualLocalPsiNo by MPI_Allgather.
657 * \param *P Problem at hand
658 * \param *Psi wave functions structure Psis
659 */
660void UpdateGramSchActualPsiNo(struct Problem *P, struct Psis *Psi)
661{
662 struct RunStruct *R = &P->R;
663 MPI_Allgather ( &(R->ActualLocalPsiNo), 1, MPI_INT,
664 Psi->AllActualLocalPsiNo, 1, MPI_INT, P->Par.comm_ST_PsiT );
665}
666
667/** Updates array Psis::AllPsiStatus from the Psis::LocalPsiStatus by MPI_Allgather.
668 * First, calculates number of MPI_OnePsiElement to be received and their displacements in the
669 * global Array Psis::AllPsiStatus, then follows MPI_Allgather and afterwards a printout to screen
670 * if verbose is specified.
671 * \param *P Problem at hand
672 * \param *Psi wave functions structure Psis
673 * \warning Don't use before FirstInitGramSch() due to needed declaration of MPI_OnePsiElement
674 */
675void UpdateGramSchAllPsiStatus(struct Problem *P, struct Psis *Psi)
676{
677 int *recvcounts, *displs;
678 int GramSchLocalNo = Psi->LocalNo+1;
679 int MyStartNo = 0;
680 int i;
681
682 //recvcounts = (int *)Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"UpdateGramSchAllPsiStatus: recvcounts");
683 //displs = (int *)Malloc(sizeof(int)*P->Par.Max_me_comm_ST_PsiT,"UpdateGramSchAllPsiStatus: displs");
684 recvcounts = (int *)Malloc(sizeof(int)*P->Par.procs,"UpdateGramSchAllPsiStatus: recvcounts");
685 displs = (int *)Malloc(sizeof(int)*P->Par.procs,"UpdateGramSchAllPsiStatus: displs");
686
687 for (i=0;i<P->Par.Max_me_comm_ST_PsiT;i++) {
688 recvcounts[i] = Psi->AllLocalNo[i]; // how many Psis should be received
689 displs[i] = MyStartNo; // displacement for these Psis
690 MyStartNo += Psi->AllLocalNo[i]; //
691 }
692 // send all (GramSchLocalNo) own local Psis and gather the AllPsiStatuses of all other processes
693 MPI_Allgatherv ( Psi->LocalPsiStatus, GramSchLocalNo, MPI_OnePsiElement,
694 Psi->AllPsiStatus, recvcounts, displs, MPI_OnePsiElement, P->Par.comm_ST_PsiT );
695
696 //if(P->Call.out[PsiOut])
697 //for (i=0;i< MyStartNo;i++)
698 //fprintf(stderr,"(%i) MyLocalNo = %i, MyGlobalNo = %i/%i, f = %lg, Type: %i, GramSch: %i, me_comm: %d, my_color_comm: %d \n",P->Par.me, Psi->AllPsiStatus[i].MyLocalNo, i, Psi->AllPsiStatus[i].MyGlobalNo, Psi->AllPsiStatus[i].PsiFactor, Psi->AllPsiStatus[i].PsiType, Psi->AllPsiStatus[i].PsiGramSchStatus, Psi->AllPsiStatus[i].me_comm_ST_Psi, Psi->AllPsiStatus[i].my_color_comm_ST_Psi);
699
700 Free(recvcounts, "UpdateGramSchAllPsiStatus: recvcounts");
701 Free(displs, "UpdateGramSchAllPsiStatus: displs");
702}
703
704/** Updates array Psis::AllOldActualLocalPsiNo from the RunStruct::OldActualLocalPsiNo by MPI_Allgather.
705 * \param *P Problem at hand
706 * \param *Psi wave functions structure Psis
707 */
708void UpdateGramSchOldActualPsiNo(struct Problem *P, struct Psis *Psi)
709{
710 struct RunStruct *R = &P->R;
711 MPI_Allgather ( &(R->OldActualLocalPsiNo), 1, MPI_INT,
712 Psi->AllOldActualLocalPsiNo, 1, MPI_INT, P->Par.comm_ST_PsiT );
713}
714
715#define max_GramSch_iter 3
716
717/** Test Gram-Schmidt-Orthogonalization.
718 * Test if all pairs of Psis are orthogonal respectively normalized (scalar product <= 1).
719 * Give output to stderr if not so.
720 * \param *P Problem at hand
721 * \param *Lev LatticeLevel structure
722 * \param *Psi wave functions structure Psis
723 * \param Type2test basically current minimisation type, see RunStruct#CurrentMin
724 */
725void TestGramSch(struct Problem *P, struct LatticeLevel *Lev, struct Psis *Psi, int Type2test) {
726 double LocalSP=0.0,PsiSP;
727 int i,j,k,s,RecvSource;
728 MPI_Status status;
729 struct OnePsiElement *OnePsiA, *LOnePsiA, *LOnePsiB;
730 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
731 int NotOk; // counts pairs that are not orthogonal
732 int iter = 0;
733 fftw_complex *LPsiDatA, *LPsiDatB;
734
735 do {
736 NotOk = 0;
737 //fprintf(stderr,"(%i) Testing Orthogonality ... \n", P->Par.me);
738 for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions (plus the extra one for each process)
739 OnePsiA = &Psi->AllPsiStatus[i]; // grab OnePsiA
740 //fprintf(stderr,"(%i) OnePsiA: Type %d, GlobalNo %d \n", P->Par.me, OnePsiA->PsiType, OnePsiA->MyGlobalNo);
741 if (OnePsiA->PsiGramSchStatus == (int)IsOrthonormal ||
742 OnePsiA->PsiGramSchStatus == (int)IsOrthogonal) { // if it has been orthogonalized
743 //fprintf(stderr,"(%i) ... orthogonal\n", P->Par.me);
744 if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
745 LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
746 else
747 LOnePsiA = NULL;
748
749 if (LOnePsiA == NULL) { // if it's not local ... receive it from respective process into TempPsi
750 RecvSource = OnePsiA->my_color_comm_ST_Psi;
751 MPI_Recv( Lev->LPsi->TempPsi, Lev->MaxG*ElementSize, MPI_DOUBLE, RecvSource, GramSchTag3, P->Par.comm_ST_PsiT, &status );
752 LPsiDatA=Lev->LPsi->TempPsi;
753 } else { // .. otherwise send it to all other processes (Max_me... - 1)
754 for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
755 if (k != OnePsiA->my_color_comm_ST_Psi)
756 MPI_Send( Lev->LPsi->LocalPsi[OnePsiA->MyLocalNo], Lev->MaxG*ElementSize, MPI_DOUBLE, k, GramSchTag3, P->Par.comm_ST_PsiT);
757 LPsiDatA=Lev->LPsi->LocalPsi[OnePsiA->MyLocalNo];
758 } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received
759
760 for (j=Psi->TypeStartIndex[Occupied]; j < Psi->TypeStartIndex[Extra]+1; j++) { // for all locally accessible including extra Psis
761 LOnePsiB = &Psi->LocalPsiStatus[j];
762 //if (LOnePsiB->PsiType > UnOccupied || OnePsiA->PsiType > UnOccupied) fprintf(stderr,"(%i) Checking global %i against local %i/%i\n",P->Par.me, OnePsiA->MyGlobalNo, LOnePsiB->MyLocalNo, LOnePsiB->MyGlobalNo);
763 if (LOnePsiB->PsiGramSchStatus == (int)IsOrthonormal ||
764 LOnePsiB->PsiGramSchStatus == (int)IsOrthogonal) { // if it's orthogonalized
765 LPsiDatB=Lev->LPsi->LocalPsi[LOnePsiB->MyLocalNo]; // set LPsiDatB onto it
766 s=0;
767 LocalSP = 0.0;
768 if (Lev->GArray[0].GSq == 0.0) { // calculate scalar product of LPsiDatA and LPsiDatB
769 LocalSP += LPsiDatA[0].re*LPsiDatB[0].re;
770 s++;
771 }
772 for (k=s; k < Lev->MaxG; k++) {
773 LocalSP += 2*(LPsiDatA[k].re*LPsiDatB[k].re+LPsiDatA[k].im*LPsiDatB[k].im);
774 }
775 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi); // gather by summation results from the group sharing the coefficients
776 //if (P->Call.out[LeaderOut])
777 switch (Type2test) {
778 default:
779 case -1: // test all, checked!
780 if (((LOnePsiB->PsiType <= UnOccupied || (LOnePsiB->MyLocalNo == P->R.ActualLocalPsiNo && OnePsiA->PsiType == Extra)) || (LOnePsiB->MyGlobalNo == OnePsiA->MyGlobalNo))) { // check if it's zero (orthogonal) and give output if wanted
781 if (i == LOnePsiB->MyGlobalNo && LOnePsiB->PsiGramSchStatus == (int)IsOrthonormal) {
782 if (fabs(PsiSP -1.0) >= MYEPSILON) {
783 fprintf(stderr,"(%i)(%i,%i) = %g ?= 1.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP -1.0), MYEPSILON);
784 NotOk++;
785 }
786 } else {
787 if (fabs(PsiSP) >= MYEPSILON && (LOnePsiB != OnePsiA && LOnePsiB->PsiType > UnOccupied)) {
788 fprintf(stderr,"(%i)(%i,%i) = %g ?= 0.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP), MYEPSILON);
789 NotOk++;
790 }
791 }
792 }
793 break;
794 case Occupied: // test unperturbed orthogonality, checked!
795 case UnOccupied:
796 if ((LOnePsiB->PsiType <= UnOccupied) && (OnePsiA->PsiType <= UnOccupied || OnePsiA->PsiType == Extra)) {
797 if (i == LOnePsiB->MyGlobalNo && LOnePsiB->PsiGramSchStatus == (int)IsOrthonormal) {
798 if (fabs(PsiSP -1.0) >= MYEPSILON) {
799 fprintf(stderr,"(%i)(%i,%i) = %g != 1.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,i,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP -1.0), MYEPSILON);
800 NotOk++;
801 } else {
802 //fprintf(stderr,"(%i)(%i,%i) = %g == 1.0 eps(%g < %g)\n",P->Par.my_color_comm_ST,i,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP -1.0), MYEPSILON);
803 }
804 } else {
805 if (fabs(PsiSP) >= MYEPSILON) {
806 fprintf(stderr,"(%i)(%i,%i) = %g != 0.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,i,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP), MYEPSILON);
807 NotOk++;
808 } else {
809 //fprintf(stderr,"(%i)(%i,%i) = %g == 0.0 eps(%g < %g)\n",P->Par.my_color_comm_ST,i,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP), MYEPSILON);
810 }
811 }
812 } else {
813 //fprintf(stderr,"(%i)(%i,%i) Not (Un)Occupied\n",P->Par.my_color_comm_ST,i,LOnePsiB->MyGlobalNo);
814 }
815 break;
816 case Perturbed_P0: // test perturbed orthogonality and normalization of all, checked!
817 case Perturbed_P1:
818 case Perturbed_P2:
819 case Perturbed_RxP0:
820 case Perturbed_RxP1:
821 case Perturbed_RxP2:
822 if ((((LOnePsiB->PsiType <= UnOccupied || LOnePsiB->PsiType == Type2test) && (OnePsiA->PsiType <= UnOccupied || OnePsiA->PsiType == Type2test) && (OnePsiA->PsiType != LOnePsiB->PsiType)) || (LOnePsiB->MyGlobalNo == OnePsiA->MyGlobalNo))) { // check if it's zero (orthogonal) and give output if wanted
823 if (i == LOnePsiB->MyGlobalNo && LOnePsiB->PsiGramSchStatus == (int)IsOrthonormal) {
824 if (fabs(PsiSP -1.0) >= MYEPSILON) {
825 fprintf(stderr,"(%i)(%i,%i) = %g != 1.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP -1.0), MYEPSILON);
826 NotOk++;
827 }
828 } else {
829 if (fabs(PsiSP) >= MYEPSILON && (LOnePsiB->MyGlobalNo != OnePsiA->MyGlobalNo)) {
830 fprintf(stderr,"(%i)(%i,%i) = %g != 0.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP), MYEPSILON);
831 NotOk++;
832 }
833 }
834 }
835 break;
836 case Extra:
837 if (((LOnePsiB->PsiType == Extra) || (LOnePsiB->PsiType == Occupied)) && ((OnePsiA->PsiType == Extra) || (OnePsiA->PsiType == Occupied))) {
838 if (i == LOnePsiB->MyGlobalNo && LOnePsiB->PsiGramSchStatus == (int)IsOrthonormal) {
839 if (fabs(PsiSP -1.0) >= MYEPSILON) {
840 fprintf(stderr,"(%i)(%i,%i) = %g != 1.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP -1.0), MYEPSILON);
841 NotOk++;
842 }
843 } else {
844 if (fabs(PsiSP) >= MYEPSILON) {
845 fprintf(stderr,"(%i)(%i,%i) = %g != 0.0 eps(%g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,LOnePsiB->MyGlobalNo,PsiSP, fabs(PsiSP), MYEPSILON);
846 NotOk++;
847 }
848 }
849 }
850 break;
851 }
852 }
853 }
854 }
855 }
856/* if (NotOk != 0) {
857 fprintf(stderr,"(%i) NotOk %i ... re-orthogonalizing type %i for the %ith time\n",P->Par.me, NotOk, Type2test, ++iter);
858 ResetGramSchTagType(P, Psi, Type2test, NotOrthogonal);
859 GramSch(P,Lev,Psi,Orthonormalize);
860 }*/
861 } while ((NotOk != 0) && (iter < max_GramSch_iter));
862 if (P->Call.out[StepLeaderOut]) { // final check if there have been any un-orthogonal pairs
863 if (Type2test != -1) {
864 if (NotOk == 0) {
865 fprintf(stderr,"(%i)TestGramSchm on %s: Ok !\n",P->Par.my_color_comm_ST, P->R.MinimisationName[Type2test]);
866 } else {
867 fprintf(stderr,"(%i)TestGramSchm on %s: There are %i pairs not Ok!\n",P->Par.my_color_comm_ST, P->R.MinimisationName[Type2test], NotOk);
868 //Error(SomeError,"Wave functions not orthonormal as they should be!");
869 }
870 } else {
871 if (NotOk == 0) {
872 fprintf(stderr,"(%i)TestGramSchm on all: Ok !\n",P->Par.my_color_comm_ST);
873 } else {
874 fprintf(stderr,"(%i)TestGramSchm on all: There are %i pairs not Ok!\n",P->Par.my_color_comm_ST,NotOk);
875 //Error(SomeError,"Wave functions not orthonormal as they should be!");
876 }
877 }
878 }
879}
880
881
882/** Test if a given wave function to all others.
883 * \param *P Problem at hand
884 * \param *Lev LatticeLevel structure
885 * \param *psi pointer to array with wave function coefficients
886 */
887void TestForOrth(struct Problem *P, struct LatticeLevel *Lev, fftw_complex *psi) {
888 struct Lattice *Lat = &P->Lat;
889 struct Psis *Psi = &Lat->Psi;
890 double LocalSP=0.0,PsiSP;
891 int i,k,RecvSource;
892 MPI_Status status;
893 struct OnePsiElement *OnePsiA, *LOnePsiA;
894 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
895 int NotOk = 0; // counts pairs that are not orthogonal
896 fftw_complex *LPsiDatA;
897 for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions (plus the extra one for each process)
898 OnePsiA = &Psi->AllPsiStatus[i]; // grab OnePsiA
899 if (OnePsiA->PsiGramSchStatus == (int)IsOrthonormal ||
900 OnePsiA->PsiGramSchStatus == (int)IsOrthogonal) { // if it has been orthogonalized
901 if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
902 LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
903 else
904 LOnePsiA = NULL;
905 if (LOnePsiA == NULL) { // if it's not local ... receive it from respective process into TempPsi
906 RecvSource = OnePsiA->my_color_comm_ST_Psi;
907 MPI_Recv( Lev->LPsi->TempPsi, Lev->MaxG*ElementSize, MPI_DOUBLE, RecvSource, GramSchTag3, P->Par.comm_ST_PsiT, &status );
908 LPsiDatA=Lev->LPsi->TempPsi;
909 } else { // .. otherwise send it to all other processes (Max_me... - 1)
910 for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
911 if (k != OnePsiA->my_color_comm_ST_Psi)
912 MPI_Send( Lev->LPsi->LocalPsi[OnePsiA->MyLocalNo], Lev->MaxG*ElementSize, MPI_DOUBLE, k, GramSchTag3, P->Par.comm_ST_PsiT);
913 LPsiDatA=Lev->LPsi->LocalPsi[OnePsiA->MyLocalNo];
914 } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received
915
916 LocalSP = 0.0;
917 k=0;
918 if (Lev->GArray[0].GSq == 0.0) { // calculate scalar product of LPsiDatA and LPsiDatB
919 LocalSP += LPsiDatA[0].re*psi[0].re;
920 k++;
921 }
922 for (; k < Lev->MaxG; k++) {
923 LocalSP += 2*(LPsiDatA[k].re*psi[k].re+LPsiDatA[k].im*psi[k].im);
924 }
925 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi); // gather by summation results from the group sharing the coefficients
926 if ((fabs(PsiSP -1.0) >= MYEPSILON) && (fabs(PsiSP) >= MYEPSILON)) {
927 fprintf(stderr,"(%i)(%i,Psi) = %g ?= 1.0 or 0.0 eps(%g or %g >= %g)\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo,PsiSP, fabs(PsiSP -1.0), fabs(PsiSP), MYEPSILON);
928 NotOk++;
929 } else
930 fprintf(stderr,"(%i)(%i,Psi) ok.\n",P->Par.my_color_comm_ST,OnePsiA->MyGlobalNo);
931 }
932 }
933 if (P->Call.out[LeaderOut]) { // final check if there have been any un-orthogonal pairs
934 if (NotOk == 0) {
935 fprintf(stderr,"(%i)TestGramSchm: Ok !\n",P->Par.my_color_comm_ST);
936 } else {
937 fprintf(stderr,"(%i)TestGramSchm: There are %i pairs not orthogonal!\n",P->Par.my_color_comm_ST,NotOk);
938 }
939 }
940}
Note: See TracBrowser for help on using the repository browser.