← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric

« back to all changes in this revision

Viewing changes to src/lib/math/scmat/disthql.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2005-11-27 11:41:49 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051127114149-zgz9r3gk50w8ww2q
Tags: 2.3.0-1
http://bugs.debian.org/339232
* New upstream release.
* debian/rules (SONAME): Activate awk snippet for automatic so-name
  detection again, resulting in a bump to `7' and making a `c2a' for
  the C++ allocator change unnecessary; closes: #339232.
* debian/patches/00list (08_gcc-4.0_fixes): Removed, no longer needed.
* debian/rules (test): Remove workarounds, do not abort build if tests
  fail.
* debian/ref: Removed.
* debian/control.in (libsc): Added Conflict against libsc6c2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/math/scmat/disthql.cc
6
6
#include <util/group/message.h>
7
7
 
8
8
#include <math/scmat/disthql.h>
9
 
#include <math/linpackd/linpackd.h>
10
9
 
11
10
#include <math/scmat/f77sym.h>
12
11
 
17
16
  void F77_PDSTEQR(int *n, double *d, double *e,
18
17
                double *z, int *ldz, int *nz, double *work,
19
18
                int *info);
 
19
  void F77_DCOPY(int *n, double *dx, int *incx, double *dy, int *incy);
 
20
  double F77_DNRM2(int *n, double *dx, int *incx);
 
21
  double F77_DDOT(int *n, double *dx, int *incx, double *dy, int *incy);
 
22
  void F77_DSCAL(int *n, double *da, double *dx, int *incx);
 
23
  void F77_DAXPY(int *n, double *da, double *dx, int *incx, double *dy, int *incy);
20
24
}
21
25
 
22
26
namespace sc {
24
28
static void dist_diagonalize_(int n, int m, double *a, double *d, double *e,
25
29
                              double *sigma, double *z, double *v, double *w,
26
30
                              int *ind, const Ref<MessageGrp>&);
27
 
static void
28
 
pimtql2_ (double *d,double *e,int *sn,double *z,int *sm,int *info);
29
 
 
30
 
static double absol(double x);
31
31
 
32
32
static void pflip(int id,int n,int m,int p,double *ar,double *ac,double *w,
33
33
                  const Ref<MessageGrp>&);
34
34
 
35
 
static double epslon (double x);
36
 
 
37
 
static void update(double *z,int m,double c,double s);
38
 
 
39
35
static void
40
36
ptred2_(double *a, int *lda, int *n, int *m, int *p, int *id,
41
37
        double *d, double *e, double *z, double *work,
108
104
 
109
105
 /* diagonalize tridiagonal matrix using implicit QL method */
110
106
 
111
 
#if 0
112
 
  pimtql2_(d,e,&n,z,&m,&info);
113
 
  if (info != 0) {
114
 
      ExEnv::outn() << "dist_diagonalize: node "
115
 
           << grp->me() << ": nonzero ierr from pimtql2" << endl;
116
 
      abort();
117
 
    }
118
 
#else
119
107
   for (i=1; i<n; i++) e[i-1] = e[i];
120
108
   F77_PDSTEQR(&n, d, e, z, &m, &m, w, &info);
121
 
#endif
122
109
 
123
110
 /* rearrange the eigenvectors by transposition */
124
111
 
125
112
  i = m * n;
126
 
  dcopy_(&i,&z[0],&one,&a[0],&one);
 
113
  F77_DCOPY(&i,&z[0],&one,&a[0],&one);
127
114
  pflip(id,n,m,nproc,&a[0],&v[0],&w[0],grp);
128
115
}
129
116
 
130
117
/* ******************************************************** */
131
 
/* Function of this subroutine :                            */ 
132
 
/*    diagonalize a real, symmetric tridiagonal matrix      */
133
 
/*    using the QL method                                   */ 
134
 
/* Parameters :                                             */
135
 
/*  on entry :                                              */
136
 
/*    d[n]    - the diagonal of the tridiagonal result      */
137
 
/*    e[n]    - the offdiagonal of the result(e[1]-e[n-1])  */
138
 
/*    sn      - size of the tridiagonal matrix              */
139
 
/*    z[m][n] - m rows of transformation matrix from before */
140
 
/*    m   - number of locally held columns                  */
141
 
/*  on return :                                             */
142
 
/*    d[n]    - the eigenvalues                             */
143
 
/*    e[n]    - non-useful information                      */
144
 
/*    z[m][n] - m rows of eigenvectors                      */
145
 
/*    info    - if 0, results are accurate                  */
146
 
/*              if nonzero, results may be inaccurate       */
147
 
/* -------------------------------------------------------- */
148
 
 
149
 
static void
150
 
pimtql2_ (double *d,double *e,int *sn,double *z,int *sm,int *info)
151
 
{
152
 
   double  c,s,t,q,u,p,h,macheps;
153
 
   int     n,m,i,j,k,im,its,maxit=30,one=1;
154
 
 
155
 
   /* extract parameters */
156
 
 
157
 
   *info = 0;
158
 
   n = *sn;
159
 
   m = *sm;
160
 
   t = 1.0;
161
 
   macheps = epslon(t);
162
 
 
163
 
   for (i=1; i<n; i++) e[i-1] = e[i];
164
 
   e[n-1] = 0.0;
165
 
   k = n - 2;
166
 
   for (j=0; j<n; j++) {
167
 
      its = 0;
168
 
      while (its < maxit) {
169
 
         for (im=j; im<=k; im++) {
170
 
            // this is the original threshold
171
 
            double threshold = macheps*(absol(d[im])+absol(d[im+1]));
172
 
            // new threshold will hopefully ensure convergence
173
 
            //if (threshold < macheps) threshold = macheps;
174
 
            if (absol(e[im]) <= threshold) break; 
175
 
//from NR:
176
 
//             double dsum = absol(d[im])+absol(d[im+1]);
177
 
//             if (dsum + absol(e[im]) == dsum) break;
178
 
         }
179
 
         u = d[j];
180
 
         if (im == j) break;
181
 
         else {
182
 
            its++;
183
 
 
184
 
            /* form implicit Wilkinson shift */
185
 
 
186
 
            q = (d[j+1] - u) / (2.0 * e[j]);
187
 
            t = sqrt(1.0 + q * q);      
188
 
            q = d[im] - u + e[j]/((q < 0.0) ? q - t : q + t);
189
 
            u = 0.0;
190
 
            s = c = 1.0;
191
 
            for (i=im-1; i>=j; i--) {
192
 
               p = s * e[i];
193
 
               h = c * e[i];
194
 
               if (absol(p) >= absol(q)) {
195
 
                  c = q / p;
196
 
                  t = sqrt(1.0 + c * c);
197
 
                  e[i+1] = p * t;
198
 
                  s = 1.0 / t;
199
 
                  c *= s;
200
 
               } else {
201
 
                  s = p / q;
202
 
                  t = sqrt(1.0 + s * s);
203
 
                  e[i+1] = q * t;
204
 
                  c = 1.0 / t;
205
 
                  s *= c;
206
 
               }
207
 
               q = d[i+1] - u;
208
 
               t = (d[i] - q) * s + 2.0 * c * h;
209
 
               u = s * t;
210
 
               d[i+1] = q + u;
211
 
               q = c * t - h;
212
 
 
213
 
               /* form eigenvectors */
214
 
 
215
 
#if 0
216
 
               for (int ia=0; ia<m; ia++) {
217
 
                  p = z[(i+1)*m+ia];
218
 
                  z[(i+1)*m+ia] = s * z[i*m+ia] + c * p; 
219
 
                  z[i*m+ia] = c * z[i*m+ia] - s * p;
220
 
               }
221
 
#else
222
 
               update(&z[i*m],m,c,s);
223
 
#endif
224
 
            }
225
 
            d[j] -= u;
226
 
            e[j] = q;
227
 
            e[im] = 0.0;
228
 
         }
229
 
      }
230
 
      if (its == maxit) {
231
 
         *info = its;
232
 
         break;
233
 
      }
234
 
   }
235
 
 
236
 
   /* order eigenvalues and eigenvectors */
237
 
 
238
 
   for (j=0; j<n-1; j++) {
239
 
      k = j;
240
 
      for (i=j+1; i<n; i++) if (d[i] < d[k]) k = i;
241
 
      if (k != j) {
242
 
         dswap_(&one,&d[j],&one,&d[k],&one);
243
 
         dswap_(&m,&z[j*m],&one,&z[k*m],&one); 
244
 
      }
245
 
   }
246
 
}
247
 
 
248
 
/* ******************************************************** */
249
 
 
250
 
static double 
251
 
absol(double x)
252
 
{
253
 
 if (x > 0.0)
254
 
   return(x);
255
 
 else
256
 
   return(-x);
257
 
}
258
 
 
259
 
/* ******************************************************** */
260
118
/* Function : transpose matrix                              */
261
119
/* -------------------------------------------------------- */
262
120
 
270
128
  for (k=0; k<n; k++) {
271
129
    r = k % p;
272
130
    if (id == r) {
273
 
      dcopy_(&n,&ar[i],&m,&w[0],&one);
 
131
      F77_DCOPY(&n,&ar[i],&m,&w[0],&one);
274
132
      i++;
275
133
    }
276
134
    grp->raw_bcast(&w[0], n*dpsize, r);
277
 
    dcopy_(&m,&w[id],&p,&ac[k],&n);
278
 
  }
279
 
}
280
 
 
281
 
/* ******************************************************** */
282
 
/* Function : calculate machine epslon                      */
283
 
/* -------------------------------------------------------- */
284
 
 
285
 
static double
286
 
epslon (double x) 
287
 
{
288
 
  double a,b,c,eps; 
289
 
 
290
 
  a = 4.0/3.0;
291
 
  eps = 0.0;
292
 
  while (eps == 0.0) {
293
 
    b = a - 1.0; 
294
 
    c = 3.0 * b; 
295
 
    eps = c-1.0; 
296
 
    if (eps < 0.0) eps = - eps;
297
 
  }
298
 
  if (x < 0.0) a = - x;
299
 
  return(eps*a); 
300
 
}
301
 
 
302
 
static void
303
 
update(double *z,int m,double c,double s)
304
 
{
305
 
  register int i;
306
 
  register double p;
307
 
 
308
 
  for (i=0; i < m; i++) {
309
 
    p = z[i+m];
310
 
    z[m+i] = s * z[i] + c * p;
311
 
    z[i]   = c * z[i] - s * p;
 
135
    F77_DCOPY(&m,&w[id],&p,&ac[k],&n);
312
136
  }
313
137
}
314
138
 
365
189
   i = sn * sm;
366
190
   alpha2 = 0.0;
367
191
   j = 0;
368
 
   dcopy_(&i,&alpha2,&j,&z[0],&inc);
 
192
   F77_DCOPY(&i,&alpha2,&j,&z[0],&inc);
369
193
   ld = sid;
370
194
   for (i=0; i<sm; i++) {
371
195
      z[ld*sm+i] = 1.0;
393
217
 
394
218
      if (sid == r) {
395
219
         q = sn - k;      
396
 
         alpha = dnrm2_(&q,&a[l*slda+k],&inc);
 
220
         alpha = F77_DNRM2(&q,&a[l*slda+k],&inc);
397
221
         if (a[l*slda+k] < 0.0) alpha = -alpha;
398
222
         if (alpha != 0.0) {
399
223
            alpha2 = 1.0 / alpha;
400
 
            dscal_(&q,&alpha2,&a[l*slda+k],&inc);
 
224
            F77_DSCAL(&q,&alpha2,&a[l*slda+k],&inc);
401
225
            a[l*slda+k] += 1.0;
402
226
         }
403
 
         dcopy_(&q,&a[l*slda+k],&inc,&d[k],&inc);
 
227
         F77_DCOPY(&q,&a[l*slda+k],&inc,&d[k],&inc);
404
228
         l++;
405
229
         ld += sp;
406
230
      }
414
238
         alpha2 = 0.0;
415
239
         j = 0;
416
240
         q = sn - k;
417
 
         dcopy_(&q,&alpha2,&j,&e[k],&inc);
 
241
         F77_DCOPY(&q,&alpha2,&j,&e[k],&inc);
418
242
         i = ld;
419
243
         for (j=l; j<sm; j++) {
420
244
            q = sn - i;
421
 
            e[i] = e[i] + ddot_(&q,&a[j*slda+i],&inc,&d[i],&inc);
 
245
            e[i] = e[i] + F77_DDOT(&q,&a[j*slda+i],&inc,&d[i],&inc);
422
246
            q = sn - i - 1;
423
 
            daxpy_(&q,&d[i],&a[slda*j+i+1],&inc,&e[i+1],&inc);
 
247
            F77_DAXPY(&q,&d[i],&a[slda*j+i+1],&inc,&e[i+1],&inc);
424
248
            i += sp;
425
249
         }
426
250
 
431
255
         if (sid == r) {
432
256
            q = sn - k;
433
257
            alpha2 = 1.0 / beta;
434
 
            dscal_(&q,&alpha2,&e[k],&inc);
435
 
            gamma = 0.5*ddot_(&q,&d[k],&inc,&e[k],&inc)/beta;
436
 
            daxpy_(&q,&gamma,&d[k],&inc,&e[k],&inc);
 
258
            F77_DSCAL(&q,&alpha2,&e[k],&inc);
 
259
            gamma = 0.5*F77_DDOT(&q,&d[k],&inc,&e[k],&inc)/beta;
 
260
            F77_DAXPY(&q,&gamma,&d[k],&inc,&e[k],&inc);
437
261
         }
438
262
 
439
263
         /* Rank two update of A, compute only lower half. */
442
266
         i = ld;
443
267
         for (j=l; j<sm; j++) {
444
268
            q = sn - i;
445
 
            daxpy_(&q,&d[i],&e[i],&inc,&a[j*slda+i],&inc);
446
 
            daxpy_(&q,&e[i],&d[i],&inc,&a[j*slda+i],&inc);
 
269
            F77_DAXPY(&q,&d[i],&e[i],&inc,&a[j*slda+i],&inc);
 
270
            F77_DAXPY(&q,&e[i],&d[i],&inc,&a[j*slda+i],&inc);
447
271
            i += sp;
448
272
         }
449
273
         q = sn - k;
450
274
         oobeta=1.0/beta;
451
275
         for (i=0; i<sm; i++) {
452
 
            gamma = ddot_(&q,&d[k],&inc,&z[k*sm+i],&sm)*oobeta;
453
 
            daxpy_(&q,&gamma,&d[k],&inc,&z[k*sm+i],&sm);
 
276
            gamma = F77_DDOT(&q,&d[k],&inc,&z[k*sm+i],&sm)*oobeta;
 
277
            F77_DAXPY(&q,&gamma,&d[k],&inc,&z[k*sm+i],&sm);
454
278
         }
455
279
      }
456
280
 
517
341
  i = sn * sm;
518
342
  alpha2 = 0.0;
519
343
  j = 0;
520
 
  dcopy_(&i, &alpha2, &j, &z[0], &inc);
 
344
  F77_DCOPY(&i, &alpha2, &j, &z[0], &inc);
521
345
  ld = sid;
522
346
  for (i = 0; i < sm; i++) {
523
347
    z[ld * sm + i] = 1.0;
549
373
      kp1l=l*slda+k;
550
374
      q = sn - k;
551
375
      atemp = a[l * slda + ld];
552
 
      alpha = dnrm2_(&q, &a[kp1l], &inc);
 
376
      alpha = F77_DNRM2(&q, &a[kp1l], &inc);
553
377
      if (a[kp1l] < 0.0) alpha = -alpha;
554
378
      if (alpha != 0.0) {
555
379
        alpha2 = 1.0 / alpha;
556
 
        dscal_(&q, &alpha2, &a[kp1l], &inc);
 
380
        F77_DSCAL(&q, &alpha2, &a[kp1l], &inc);
557
381
        a[kp1l] += 1.0;
558
382
        }
559
383
 
569
393
 
570
394
        if (beta != 0.0) {
571
395
          for (i = 0; i < sm; i++) {
572
 
            gamma = ddot_(&nmik, &d[ik], &inc, &z[ik * sm + i], &sm) / beta;
573
 
            daxpy_(&nmik, &gamma, &d[ik], &inc, &z[ik * sm + i], &sm);
 
396
            gamma = F77_DDOT(&nmik, &d[ik], &inc, &z[ik * sm + i], &sm) / beta;
 
397
            F77_DAXPY(&nmik, &gamma, &d[ik], &inc, &z[ik * sm + i], &sm);
574
398
            }
575
399
          }
576
400
        e[ik] = 0.0;
578
402
        }
579
403
 
580
404
   /* now resume normal service */
581
 
      dcopy_(&q, &a[kp1l], &inc, &d[k], &inc);
 
405
      F77_DCOPY(&q, &a[kp1l], &inc, &d[k], &inc);
582
406
      l++;
583
407
      ld += sp;
584
408
      }
595
419
      alpha2 = 0.0;
596
420
      j = 0;
597
421
      q = sn - k;
598
 
      dcopy_(&q, &alpha2, &j, &e[k], &inc);
 
422
      F77_DCOPY(&q, &alpha2, &j, &e[k], &inc);
599
423
      i = ld;
600
424
      for (j = l; j < sm; j++) {
601
425
        int ij=j*slda+i;
602
426
        q = sn - i;
603
 
        e[i] += ddot_(&q, &a[ij], &inc, &d[i], &inc);
 
427
        e[i] += F77_DDOT(&q, &a[ij], &inc, &d[i], &inc);
604
428
        q--;
605
 
        daxpy_(&q, &d[i], &a[ij+1], &inc, &e[i + 1], &inc);
 
429
        F77_DAXPY(&q, &d[i], &a[ij+1], &inc, &e[i + 1], &inc);
606
430
        i += sp;
607
431
        }
608
432
      grp->sum(&e[k], sn-k, work);
614
438
 
615
439
      q = sn - k;
616
440
      alpha2 = 1.0 / beta;
617
 
      dscal_(&q, &alpha2, &e[k], &inc);
618
 
      gamma = 0.5 * ddot_(&q, &d[k], &inc, &e[k], &inc) / beta;
619
 
      daxpy_(&q, &gamma, &d[k], &inc, &e[k], &inc);
 
441
      F77_DSCAL(&q, &alpha2, &e[k], &inc);
 
442
      gamma = 0.5 * F77_DDOT(&q, &d[k], &inc, &e[k], &inc) / beta;
 
443
      F77_DAXPY(&q, &gamma, &d[k], &inc, &e[k], &inc);
620
444
 
621
445
      /* Rank two update of A, compute only lower half. */
622
446
      /* A  =  A + u'*v + v'*u  =  H*A*H                */
625
449
      for (j = l; j < sm; j++) {
626
450
        double *atmp= &a[j*slda+i];
627
451
        q = sn - i;
628
 
        daxpy_(&q, &d[i], &e[i], &inc, atmp, &inc);
629
 
        daxpy_(&q, &e[i], &d[i], &inc, atmp, &inc);
 
452
        F77_DAXPY(&q, &d[i], &e[i], &inc, atmp, &inc);
 
453
        F77_DAXPY(&q, &e[i], &d[i], &inc, atmp, &inc);
630
454
        i += sp;
631
455
        }
632
456
 
640
464
        q = sn - k;
641
465
        oobeta = 1.0 / beta;
642
466
        for (i = 0; i < sm; i++) {
643
 
          gamma = ddot_(&q, &d[k], &inc, &z[k * sm + i], &sm) * oobeta;
644
 
          daxpy_(&q, &gamma, &d[k], &inc, &z[k * sm + i], &sm);
 
467
          gamma = F77_DDOT(&q, &d[k], &inc, &z[k * sm + i], &sm) * oobeta;
 
468
          F77_DAXPY(&q, &gamma, &d[k], &inc, &z[k * sm + i], &sm);
645
469
          }
646
470
        }
647
471
      }