Actual source code: dlaed3m.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    BDC - Block-divide and conquer (see description in README file)
 12: */

 14: #include <slepc/private/dsimpl.h>
 15: #include <slepcblaslapack.h>

 17: PetscErrorCode BDC_dlaed3m_(const char *jobz,const char *defl,PetscBLASInt k,PetscBLASInt n,
 18:         PetscBLASInt n1,PetscReal *d,PetscReal *q,PetscBLASInt ldq,
 19:         PetscReal rho,PetscReal *dlamda,PetscReal *q2,PetscBLASInt *indx,
 20:         PetscBLASInt *ctot,PetscReal *w,PetscReal *s,PetscBLASInt *info,
 21:         PetscBLASInt jobz_len,PetscBLASInt defl_len)
 22: {
 23: /*  -- Routine written in LAPACK version 3.0 style -- */
 24: /* *************************************************** */
 25: /*     Written by */
 26: /*     Michael Moldaschl and Wilfried Gansterer */
 27: /*     University of Vienna */
 28: /*     last modification: March 16, 2014 */

 30: /*     Small adaptations of original code written by */
 31: /*     Wilfried Gansterer and Bob Ward, */
 32: /*     Department of Computer Science, University of Tennessee */
 33: /*     see https://doi.org/10.1137/S1064827501399432 */
 34: /* *************************************************** */

 36: /*  Purpose */
 37: /*  ======= */

 39: /*  DLAED3M finds the roots of the secular equation, as defined by the */
 40: /*  values in D, W, and RHO, between 1 and K.  It makes the */
 41: /*  appropriate calls to DLAED4 and then updates the eigenvectors by */
 42: /*  multiplying the matrix of eigenvectors of the pair of eigensystems */
 43: /*  being combined by the matrix of eigenvectors of the K-by-K system */
 44: /*  which is solved here. */

 46: /*  This code makes very mild assumptions about floating point */
 47: /*  arithmetic. It will work on machines with a guard digit in */
 48: /*  add/subtract, or on those binary machines without guard digits */
 49: /*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. */
 50: /*  It could conceivably fail on hexadecimal or decimal machines */
 51: /*  without guard digits, but we know of none. */

 53: /*  Arguments */
 54: /*  ========= */

 56: /*  JOBZ    (input) CHARACTER*1 */
 57: /*          = 'N':  Do not accumulate eigenvectors (not implemented); */
 58: /*          = 'D':  Do accumulate eigenvectors in the divide-and-conquer */
 59: /*                  process. */

 61: /*  DEFL    (input) CHARACTER*1 */
 62: /*          = '0':  No deflation happened in DSRTDF */
 63: /*          = '1':  Some deflation happened in DSRTDF (and therefore some */
 64: /*                  Givens rotations need to be applied to the computed */
 65: /*                  eigenvector matrix Q) */

 67: /*  K       (input) INTEGER */
 68: /*          The number of terms in the rational function to be solved by */
 69: /*          DLAED4. 0 <= K <= N. */

 71: /*  N       (input) INTEGER */
 72: /*          The number of rows and columns in the Q matrix. */
 73: /*          N >= K (deflation may result in N>K). */

 75: /*  N1      (input) INTEGER */
 76: /*          The location of the last eigenvalue in the leading submatrix. */
 77: /*          min(1,N) <= N1 <= max(1,N-1). */

 79: /*  D       (output) DOUBLE PRECISION array, dimension (N) */
 80: /*          D(I) contains the updated eigenvalues for */
 81: /*          1 <= I <= K. */

 83: /*  Q       (output) DOUBLE PRECISION array, dimension (LDQ,N) */
 84: /*          Initially the first K columns are used as workspace. */
 85: /*          On output the columns 1 to K contain */
 86: /*          the updated eigenvectors. */

 88: /*  LDQ     (input) INTEGER */
 89: /*          The leading dimension of the array Q.  LDQ >= max(1,N). */

 91: /*  RHO     (input) DOUBLE PRECISION */
 92: /*          The value of the parameter in the rank one update equation. */
 93: /*          RHO >= 0 required. */

 95: /*  DLAMDA  (input/output) DOUBLE PRECISION array, dimension (K) */
 96: /*          The first K elements of this array contain the old roots */
 97: /*          of the deflated updating problem.  These are the poles */
 98: /*          of the secular equation. May be changed on output by */
 99: /*          having lowest order bit set to zero on Cray X-MP, Cray Y-MP, */
100: /*          Cray-2, or Cray C-90, as described above. */

102: /*  Q2      (input) DOUBLE PRECISION array, dimension (LDQ2, N) */
103: /*          The first K columns of this matrix contain the non-deflated */
104: /*          eigenvectors for the split problem. */

106: /*  INDX    (input) INTEGER array, dimension (N) */
107: /*          The permutation used to arrange the columns of the deflated */
108: /*          Q matrix into three groups (see DLAED2). */
109: /*          The rows of the eigenvectors found by DLAED4 must be likewise */
110: /*          permuted before the matrix multiply can take place. */

112: /*  CTOT    (input) INTEGER array, dimension (4) */
113: /*          A count of the total number of the various types of columns */
114: /*          in Q, as described in INDX.  The fourth column type is any */
115: /*          column which has been deflated. */

117: /*  W       (input/output) DOUBLE PRECISION array, dimension (K) */
118: /*          The first K elements of this array contain the components */
119: /*          of the deflation-adjusted updating vector. Destroyed on */
120: /*          output. */

122: /*  S       (workspace) DOUBLE PRECISION array, dimension */
123: /*          (MAX(CTOT(1)+CTOT(2),CTOT(2)+CTOT(3)) + 1)*K */
124: /*          Will contain parts of the eigenvectors of the repaired matrix */
125: /*          which will be multiplied by the previously accumulated */
126: /*          eigenvectors to update the system. This array is a major */
127: /*          source of workspace requirements ! */

129: /*  INFO    (output) INTEGER */
130: /*          = 0:  successful exit. */
131: /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
132: /*          > 0:  if INFO = i, eigenpair i was not computed successfully */

134: /*  Further Details */
135: /*  =============== */

137: /*  Based on code written by */
138: /*     Wilfried Gansterer and Bob Ward, */
139: /*     Department of Computer Science, University of Tennessee */
140: /*  Based on the design of the LAPACK code DLAED3 with small modifications */
141: /*  (Note that in contrast to the original DLAED3, this routine */
142: /*  DOES NOT require that N1 <= N/2) */

144: /*  Based on contributions by */
145: /*     Jeff Rutter, Computer Science Division, University of California */
146: /*     at Berkeley, USA */
147: /*  Modified by Francoise Tisseur, University of Tennessee. */

149: /*  ===================================================================== */

151:   PetscReal    temp, done = 1.0, dzero = 0.0;
152:   PetscBLASInt i, j, n2, n12, ii, n23, iq2, i1, one=1;

154:   *info = 0;

156:   if (k < 0) *info = -3;
157:   else if (n < k) *info = -4;
158:   else if (n1 < PetscMin(1,n) || n1 > PetscMax(1,n)) *info = -5;
159:   else if (ldq < PetscMax(1,n)) *info = -8;
160:   else if (rho < 0.) *info = -9;

163:   /* Quick return if possible */

165:   if (k == 0) PetscFunctionReturn(0);

167:   /* Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can */
168:   /* be computed with high relative accuracy (barring over/underflow). */
169:   /* This is a problem on machines without a guard digit in */
170:   /* add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). */
171:   /* The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), */
172:   /* which on any of these machines zeros out the bottommost */
173:   /* bit of DLAMDA(I) if it is 1; this makes the subsequent */
174:   /* subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation */
175:   /* occurs. On binary machines with a guard digit (almost all */
176:   /* machines) it does not change DLAMDA(I) at all. On hexadecimal */
177:   /* and decimal machines with a guard digit, it slightly */
178:   /* changes the bottommost bits of DLAMDA(I). It does not account */
179:   /* for hexadecimal or decimal machines without guard digits */
180:   /* (we know of none). We use a subroutine call to compute */
181:   /* 2*DLAMBDA(I) to prevent optimizing compilers from eliminating */
182:   /* this code. */

184:   for (i = 0; i < k; ++i) {
185:     dlamda[i] = LAPACKlamc3_(&dlamda[i], &dlamda[i]) - dlamda[i];
186:   }

188:   for (j = 1; j <= k; ++j) {

190:     /* ....calling DLAED4 for eigenpair J.... */

192:     PetscStackCallBLAS("LAPACKlaed4",LAPACKlaed4_(&k, &j, dlamda, w, &q[(j-1)*ldq], &rho, &d[j-1], info));
193:     SlepcCheckLapackInfo("laed4",*info);

195:     if (j < k) {

197:       /* If the zero finder terminated properly, but the computed */
198:       /* eigenvalues are not ordered, issue an error statement */
199:       /* but continue computation. */

203:     }
204:   }

206:   if (k == 1) goto L110;

208:   if (k == 2) {

210:     /* permute the components of Q(:,J) (the information returned by DLAED4 */
211:     /* necessary to construct the eigenvectors) according to the permutation */
212:     /* stored in INDX, resulting from deflation */

214:     for (j = 0; j < k; ++j) {
215:       w[0] = q[0+j*ldq];
216:       w[1] = q[1+j*ldq];
217:       ii = indx[0];
218:       q[0+j*ldq] = w[ii-1];
219:       ii = indx[1];
220:       q[1+j*ldq] = w[ii-1];
221:     }
222:     goto L110;
223:   }

225:   /* ....K.GE.3.... */
226:   /* Compute updated W (used for computing the eigenvectors corresponding */
227:   /* to the previously computed eigenvalues). */

229:   PetscStackCallBLAS("BLAScopy",BLAScopy_(&k, w, &one, s, &one));

231:   /* Initialize W(I) = Q(I,I) */

233:   i1 = ldq + 1;
234:   PetscStackCallBLAS("BLAScopy",BLAScopy_(&k, q, &i1, w, &one));
235:   for (j = 0; j < k; ++j) {
236:     for (i = 0; i < j; ++i) {
237:       w[i] *= q[i+j*ldq] / (dlamda[i] - dlamda[j]);
238:     }
239:     for (i = j + 1; i < k; ++i) {
240:       w[i] *= q[i+j*ldq] / (dlamda[i] - dlamda[j]);
241:     }
242:   }
243:   for (i = 0; i < k; ++i) {
244:     temp = PetscSqrtReal(-w[i]);
245:     if (temp<0) temp = -temp;
246:     w[i] =  (s[i] >= 0) ? temp : -temp;
247:   }

249:   /* Compute eigenvectors of the modified rank-1 modification (using the */
250:   /* vector W). */

252:   for (j = 0; j < k; ++j) {
253:     for (i = 0; i < k; ++i) {
254:       s[i] = w[i] / q[i+j*ldq];
255:     }
256:     temp = BLASnrm2_(&k, s, &one);
257:     for (i = 0; i < k; ++i) {

259:       /* apply the permutation resulting from deflation as stored */
260:       /* in INDX */

262:       ii = indx[i];
263:       q[i+j*ldq] = s[ii-1] / temp;
264:     }
265:   }

267: /* ************************************************************************** */

269:   /* ....updating the eigenvectors.... */

271: L110:

273:   n2 = n - n1;
274:   n12 = ctot[0] + ctot[1];
275:   n23 = ctot[1] + ctot[2];
276:   if (*(unsigned char *)jobz == 'D') {

278:     /* Compute the updated eigenvectors. (NOTE that every call of */
279:     /* DGEMM requires three DISTINCT arrays) */

281:     /* copy Q(CTOT(1)+1:K,1:K) to S */

283:     for (j=0;j<k;j++) for (i=0;i<n23;i++) s[i+j*n23] = q[ctot[0]+i+j*ldq];
284:     iq2 = n1 * n12 + 1;

286:     if (n23 != 0) {

288:       /* multiply the second part of Q2 (the eigenvectors of the */
289:       /* lower block) with S and write the result into the lower part of */
290:       /* Q, i.e., Q(N1+1:N,1:K) */

292:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N", "N", &n2, &k, &n23, &done,
293:                   &q2[iq2-1], &n2, s, &n23, &dzero, &q[n1], &ldq));
294:     } else {
295:       for (j=0;j<k;j++) for (i=0;i<n2;i++) q[n1+i+j*ldq] = 0.0;
296:     }

298:     /* copy Q(1:CTOT(1)+CTOT(2),1:K) to S */

300:     for (j=0;j<k;j++) for (i=0;i<n12;i++) s[i+j*n12] = q[i+j*ldq];

302:     if (n12 != 0) {

304:       /* multiply the first part of Q2 (the eigenvectors of the */
305:       /* upper block) with S and write the result into the upper part of */
306:       /* Q, i.e., Q(1:N1,1:K) */

308:       PetscStackCallBLAS("BLASgemm",BLASgemm_("N", "N", &n1, &k, &n12, &done,
309:                   q2, &n1, s, &n12, &dzero, q, &ldq));
310:     } else {
311:       for (j=0;j<k;j++) for (i=0;i<n1;i++) q[i+j*ldq] = 0.0;
312:     }
313:   }
314:   PetscFunctionReturn(0);
315: }