Actual source code: dspriv.c

slepc-3.15.0 2021-03-31
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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:    Private DS routines
 12: */

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

 17: PetscErrorCode DSAllocateMatrix_Private(DS ds,DSMatType m,PetscBool isreal)
 18: {
 19:   size_t         sz;
 20:   PetscInt       n,d,nelem;
 21:   PetscBool      ispep;

 25:   PetscObjectTypeCompare((PetscObject)ds,DSPEP,&ispep);
 26:   if (ispep) {
 27:     DSPEPGetDegree(ds,&d);
 28:   }
 29:   if (ispep && (m==DS_MAT_A || m==DS_MAT_B || m==DS_MAT_W || m==DS_MAT_U || m==DS_MAT_X || m==DS_MAT_Y)) n = d*ds->ld;
 30:   else n = ds->ld;
 31:   switch (m) {
 32:     case DS_MAT_T:
 33:       nelem = 3*ds->ld;
 34:       break;
 35:     case DS_MAT_D:
 36:       nelem = ds->ld;
 37:       break;
 38:     case DS_MAT_X:
 39:       nelem = ds->ld*n;
 40:       break;
 41:     case DS_MAT_Y:
 42:       nelem = ds->ld*n;
 43:       break;
 44:     default:
 45:       nelem = n*n;
 46:   }
 47:   if (isreal) {
 48:     sz = nelem*sizeof(PetscReal);
 49:     if (ds->rmat[m]) {
 50:       PetscFree(ds->rmat[m]);
 51:     } else {
 52:       PetscLogObjectMemory((PetscObject)ds,sz);
 53:     }
 54:     PetscCalloc1(nelem,&ds->rmat[m]);
 55:   } else {
 56:     sz = nelem*sizeof(PetscScalar);
 57:     if (ds->mat[m]) {
 58:       PetscFree(ds->mat[m]);
 59:     } else {
 60:       PetscLogObjectMemory((PetscObject)ds,sz);
 61:     }
 62:     PetscCalloc1(nelem,&ds->mat[m]);
 63:   }
 64:   return(0);
 65: }

 67: PetscErrorCode DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
 68: {

 72:   if (s>ds->lwork) {
 73:     PetscFree(ds->work);
 74:     PetscMalloc1(s,&ds->work);
 75:     PetscLogObjectMemory((PetscObject)ds,(s-ds->lwork)*sizeof(PetscScalar));
 76:     ds->lwork = s;
 77:   }
 78:   if (r>ds->lrwork) {
 79:     PetscFree(ds->rwork);
 80:     PetscMalloc1(r,&ds->rwork);
 81:     PetscLogObjectMemory((PetscObject)ds,(r-ds->lrwork)*sizeof(PetscReal));
 82:     ds->lrwork = r;
 83:   }
 84:   if (i>ds->liwork) {
 85:     PetscFree(ds->iwork);
 86:     PetscMalloc1(i,&ds->iwork);
 87:     PetscLogObjectMemory((PetscObject)ds,(i-ds->liwork)*sizeof(PetscBLASInt));
 88:     ds->liwork = i;
 89:   }
 90:   return(0);
 91: }

 93: /*@C
 94:    DSViewMat - Prints one of the internal DS matrices.

 96:    Collective on ds

 98:    Input Parameters:
 99: +  ds     - the direct solver context
100: .  viewer - visualization context
101: -  m      - matrix to display

103:    Note:
104:    Works only for ascii viewers. Set the viewer in Matlab format if
105:    want to paste into Matlab.

107:    Level: developer

109: .seealso: DSView()
110: @*/
111: PetscErrorCode DSViewMat(DS ds,PetscViewer viewer,DSMatType m)
112: {
113:   PetscErrorCode    ierr;
114:   PetscInt          i,j,rows,cols;
115:   PetscScalar       *v;
116:   PetscViewerFormat format;
117: #if defined(PETSC_USE_COMPLEX)
118:   PetscBool         allreal = PETSC_TRUE;
119: #endif

124:   DSCheckValidMat(ds,m,3);
125:   if (!viewer) {
126:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer);
127:   }
130:   PetscViewerGetFormat(viewer,&format);
131:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
132:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
133:   DSMatGetSize(ds,m,&rows,&cols);
134: #if defined(PETSC_USE_COMPLEX)
135:   /* determine if matrix has all real values */
136:   v = ds->mat[m];
137:   for (i=0;i<rows;i++)
138:     for (j=0;j<cols;j++)
139:       if (PetscImaginaryPart(v[i+j*ds->ld])) { allreal = PETSC_FALSE; break; }
140: #endif
141:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
142:     PetscViewerASCIIPrintf(viewer,"%% Size = %D %D\n",rows,cols);
143:     PetscViewerASCIIPrintf(viewer,"%s = [\n",DSMatName[m]);
144:   } else {
145:     PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[m]);
146:   }

148:   for (i=0;i<rows;i++) {
149:     v = ds->mat[m]+i;
150:     for (j=0;j<cols;j++) {
151: #if defined(PETSC_USE_COMPLEX)
152:       if (allreal) {
153:         PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
154:       } else {
155:         PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
156:       }
157: #else
158:       PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
159: #endif
160:       v += ds->ld;
161:     }
162:     PetscViewerASCIIPrintf(viewer,"\n");
163:   }

165:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
166:     PetscViewerASCIIPrintf(viewer,"];\n");
167:   }
168:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
169:   PetscViewerFlush(viewer);
170:   return(0);
171: }

173: PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
174: {
176:   PetscScalar    re,im,wi0;
177:   PetscInt       n,i,j,result,tmp1,tmp2=0,d=1;

180:   n = ds->t;   /* sort only first t pairs if truncated */
181:   /* insertion sort */
182:   i=ds->l+1;
183: #if !defined(PETSC_USE_COMPLEX)
184:   if (wi && wi[perm[i-1]]!=0.0) i++; /* initial value is complex */
185: #else
186:   if (isghiep && PetscImaginaryPart(wr[perm[i-1]])!=0.0) i++;
187: #endif
188:   for (;i<n;i+=d) {
189:     re = wr[perm[i]];
190:     if (wi) im = wi[perm[i]];
191:     else im = 0.0;
192:     tmp1 = perm[i];
193: #if !defined(PETSC_USE_COMPLEX)
194:     if (im!=0.0) { d = 2; tmp2 = perm[i+1]; }
195:     else d = 1;
196: #else
197:     if (isghiep && PetscImaginaryPart(re)!=0.0) { d = 2; tmp2 = perm[i+1]; }
198:     else d = 1;
199: #endif
200:     j = i-1;
201:     if (wi) wi0 = wi[perm[j]];
202:     else wi0 = 0.0;
203:     SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result);
204:     while (result<0 && j>=ds->l) {
205:       perm[j+d] = perm[j];
206:       j--;
207: #if !defined(PETSC_USE_COMPLEX)
208:       if (wi && wi[perm[j+1]]!=0)
209: #else
210:       if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
211: #endif
212:         { perm[j+d] = perm[j]; j--; }

214:       if (j>=ds->l) {
215:         if (wi) wi0 = wi[perm[j]];
216:         else wi0 = 0.0;
217:         SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result);
218:       }
219:     }
220:     perm[j+1] = tmp1;
221:     if (d==2) perm[j+2] = tmp2;
222:   }
223:   return(0);
224: }

226: PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
227: {
229:   PetscScalar    re;
230:   PetscInt       i,j,result,tmp,l,n;

233:   n = ds->t;   /* sort only first t pairs if truncated */
234:   l = ds->l;
235:   /* insertion sort */
236:   for (i=l+1;i<n;i++) {
237:     re = eig[perm[i]];
238:     j = i-1;
239:     SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result);
240:     while (result<0 && j>=l) {
241:       tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
242:       if (j>=l) {
243:         SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result);
244:       }
245:     }
246:   }
247:   return(0);
248: }

250: /*
251:   DSCopyMatrix_Private - Copies the trailing block of a matrix (from
252:   rows/columns l to n).
253: */
254: PetscErrorCode DSCopyMatrix_Private(DS ds,DSMatType dst,DSMatType src)
255: {
257:   PetscInt    j,m,off,ld;
258:   PetscScalar *S,*D;

261:   ld  = ds->ld;
262:   m   = ds->n-ds->l;
263:   off = ds->l+ds->l*ld;
264:   S   = ds->mat[src];
265:   D   = ds->mat[dst];
266:   for (j=0;j<m;j++) {
267:     PetscArraycpy(D+off+j*ld,S+off+j*ld,m);
268:   }
269:   return(0);
270: }

272: PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt l,PetscInt n,DSMatType mat,PetscInt *perm)
273: {
274:   PetscInt    i,j,k,p,ld;
275:   PetscScalar *Q,rtmp;

278:   ld = ds->ld;
279:   Q  = ds->mat[mat];
280:   for (i=l;i<n;i++) {
281:     p = perm[i];
282:     if (p != i) {
283:       j = i + 1;
284:       while (perm[j] != i) j++;
285:       perm[j] = p; perm[i] = i;
286:       /* swap columns i and j */
287:       for (k=0;k<n;k++) {
288:         rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
289:       }
290:     }
291:   }
292:   return(0);
293: }

295: PetscErrorCode DSPermuteRows_Private(DS ds,PetscInt l,PetscInt n,DSMatType mat,PetscInt *perm)
296: {
297:   PetscInt    i,j,m=ds->m,k,p,ld;
298:   PetscScalar *Q,rtmp;

301:   if (!m) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"m was not set");
302:   ld = ds->ld;
303:   Q  = ds->mat[mat];
304:   for (i=l;i<n;i++) {
305:     p = perm[i];
306:     if (p != i) {
307:       j = i + 1;
308:       while (perm[j] != i) j++;
309:       perm[j] = p; perm[i] = i;
310:       /* swap rows i and j */
311:       for (k=0;k<m;k++) {
312:         rtmp = Q[p+k*ld]; Q[p+k*ld] = Q[i+k*ld]; Q[i+k*ld] = rtmp;
313:       }
314:     }
315:   }
316:   return(0);
317: }

319: PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt l,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
320: {
321:   PetscInt    i,j,m=ds->m,k,p,ld;
322:   PetscScalar *U,*VT,rtmp;

325:   if (!m) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"m was not set");
326:   ld = ds->ld;
327:   U  = ds->mat[mat1];
328:   VT = ds->mat[mat2];
329:   for (i=l;i<n;i++) {
330:     p = perm[i];
331:     if (p != i) {
332:       j = i + 1;
333:       while (perm[j] != i) j++;
334:       perm[j] = p; perm[i] = i;
335:       /* swap columns i and j of U */
336:       for (k=0;k<n;k++) {
337:         rtmp = U[k+p*ld]; U[k+p*ld] = U[k+i*ld]; U[k+i*ld] = rtmp;
338:       }
339:       /* swap rows i and j of VT */
340:       for (k=0;k<m;k++) {
341:         rtmp = VT[p+k*ld]; VT[p+k*ld] = VT[i+k*ld]; VT[i+k*ld] = rtmp;
342:       }
343:     }
344:   }
345:   return(0);
346: }

348: /*@
349:    DSSetIdentity - Copy the identity (a diagonal matrix with ones) on the
350:    active part of a matrix.

352:    Logically Collective on ds

354:    Input Parameters:
355: +  ds  - the direct solver context
356: -  mat - the matrix to modify

358:    Level: intermediate
359: @*/
360: PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
361: {
363:   PetscScalar    *x;
364:   PetscInt       i,ld,n,l;

369:   DSCheckValidMat(ds,mat,2);

371:   DSGetDimensions(ds,&n,NULL,&l,NULL,NULL);
372:   DSGetLeadingDimension(ds,&ld);
373:   PetscLogEventBegin(DS_Other,ds,0,0,0);
374:   DSGetArray(ds,mat,&x);
375:   PetscArrayzero(&x[ld*l],ld*(n-l));
376:   for (i=l;i<n;i++) x[i+i*ld] = 1.0;
377:   DSRestoreArray(ds,mat,&x);
378:   PetscLogEventEnd(DS_Other,ds,0,0,0);
379:   return(0);
380: }

382: /*@C
383:    DSOrthogonalize - Orthogonalize the columns of a matrix.

385:    Logically Collective on ds

387:    Input Parameters:
388: +  ds   - the direct solver context
389: .  mat  - a matrix
390: -  cols - number of columns to orthogonalize (starting from column zero)

392:    Output Parameter:
393: .  lindcols - (optional) number of linearly independent columns of the matrix

395:    Level: developer

397: .seealso: DSPseudoOrthogonalize()
398: @*/
399: PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
400: {
402:   PetscInt       n,l,ld;
403:   PetscBLASInt   ld_,rA,cA,info,ltau,lw;
404:   PetscScalar    *A,*tau,*w,saux,dummy;

408:   DSCheckAlloc(ds,1);
410:   DSCheckValidMat(ds,mat,2);

413:   DSGetDimensions(ds,&n,NULL,&l,NULL,NULL);
414:   DSGetLeadingDimension(ds,&ld);
415:   n = n - l;
416:   if (cols > n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
417:   if (n == 0 || cols == 0) return(0);

419:   PetscLogEventBegin(DS_Other,ds,0,0,0);
420:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
421:   DSGetArray(ds,mat,&A);
422:   PetscBLASIntCast(PetscMin(cols,n),&ltau);
423:   PetscBLASIntCast(ld,&ld_);
424:   PetscBLASIntCast(n,&rA);
425:   PetscBLASIntCast(cols,&cA);
426:   lw = -1;
427:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
428:   SlepcCheckLapackInfo("geqrf",info);
429:   lw = (PetscBLASInt)PetscRealPart(saux);
430:   DSAllocateWork_Private(ds,lw+ltau,0,0);
431:   tau = ds->work;
432:   w = &tau[ltau];
433:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
434:   SlepcCheckLapackInfo("geqrf",info);
435:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,&ltau,&ltau,&A[ld*l+l],&ld_,tau,w,&lw,&info));
436:   SlepcCheckLapackInfo("orgqr",info);
437:   if (lindcols) *lindcols = ltau;

439:   PetscFPTrapPop();
440:   PetscLogEventEnd(DS_Other,ds,0,0,0);
441:   DSRestoreArray(ds,mat,&A);
442:   PetscObjectStateIncrease((PetscObject)ds);
443:   return(0);
444: }

446: /*
447:   Compute C <- a*A*B + b*C, where
448:     ldC, the leading dimension of C,
449:     ldA, the leading dimension of A,
450:     rA, cA, rows and columns of A,
451:     At, if true use the transpose of A instead,
452:     ldB, the leading dimension of B,
453:     rB, cB, rows and columns of B,
454:     Bt, if true use the transpose of B instead
455: */
456: static PetscErrorCode SlepcMatDenseMult(PetscScalar *C,PetscInt _ldC,PetscScalar b,PetscScalar a,const PetscScalar *A,PetscInt _ldA,PetscInt rA,PetscInt cA,PetscBool At,const PetscScalar *B,PetscInt _ldB,PetscInt rB,PetscInt cB,PetscBool Bt)
457: {
459:   PetscInt       tmp;
460:   PetscBLASInt   m, n, k, ldA = _ldA, ldB = _ldB, ldC = _ldC;
461:   const char     *N = "N", *T = "C", *qA = N, *qB = N;

464:   if ((rA == 0) || (cB == 0)) return(0);

469:   /* Transpose if needed */
470:   if (At) tmp = rA, rA = cA, cA = tmp, qA = T;
471:   if (Bt) tmp = rB, rB = cB, cB = tmp, qB = T;

473:   /* Check size */
474:   if (cA != rB) SETERRQ(PETSC_COMM_SELF,1,"Matrix dimensions do not match");

476:   /* Do stub */
477:   if ((rA == 1) && (cA == 1) && (cB == 1)) {
478:     if (!At && !Bt) *C = *A * *B;
479:     else if (At && !Bt) *C = PetscConj(*A) * *B;
480:     else if (!At && Bt) *C = *A * PetscConj(*B);
481:     else *C = PetscConj(*A) * PetscConj(*B);
482:     m = n = k = 1;
483:   } else {
484:     m = rA; n = cB; k = cA;
485:     PetscStackCallBLAS("BLASgemm",BLASgemm_(qA,qB,&m,&n,&k,&a,(PetscScalar*)A,&ldA,(PetscScalar*)B,&ldB,&b,C,&ldC));
486:   }

488:   PetscLogFlops(2.0*m*n*k);
489:   return(0);
490: }

492: /*@C
493:    DSPseudoOrthogonalize - Orthogonalize the columns of a matrix with Modified
494:    Gram-Schmidt in an indefinite inner product space defined by a signature.

496:    Logically Collective on ds

498:    Input Parameters:
499: +  ds   - the direct solver context
500: .  mat  - the matrix
501: .  cols - number of columns to orthogonalize (starting from column zero)
502: -  s    - the signature that defines the inner product

504:    Output Parameters:
505: +  lindcols - (optional) linearly independent columns of the matrix
506: -  ns   - (optional) the new signature of the vectors

508:    Note:
509:    After the call the matrix satisfies A'*s*A = ns.

511:    Level: developer

513: .seealso: DSOrthogonalize()
514: @*/
515: PetscErrorCode DSPseudoOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscReal *s,PetscInt *lindcols,PetscReal *ns)
516: {
518:   PetscInt       i,j,k,l,n,ld;
519:   PetscBLASInt   info,one=1,zero=0,rA_,ld_;
520:   PetscScalar    *A,*A_,*m,*h,nr0;
521:   PetscReal      nr_o,nr,nr_abs,*ns_,done=1.0;

525:   DSCheckAlloc(ds,1);
527:   DSCheckValidMat(ds,mat,2);
530:   DSGetDimensions(ds,&n,NULL,&l,NULL,NULL);
531:   DSGetLeadingDimension(ds,&ld);
532:   n = n - l;
533:   if (cols > n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
534:   if (n == 0 || cols == 0) return(0);
535:   PetscBLASIntCast(n,&rA_);
536:   PetscBLASIntCast(ld,&ld_);
537:   DSGetArray(ds,mat,&A_);
538:   A = &A_[ld*l+l];
539:   DSAllocateWork_Private(ds,n+cols,ns?0:cols,0);
540:   m = ds->work;
541:   h = &m[n];
542:   ns_ = ns ? ns : ds->rwork;
543:   PetscLogEventBegin(DS_Other,ds,0,0,0);
544:   for (i=0; i<cols; i++) {
545:     /* m <- diag(s)*A[i] */
546:     for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
547:     /* nr_o <- mynorm(A[i]'*m), mynorm(x) = sign(x)*sqrt(|x|) */
548:     SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
549:     nr = nr_o = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
550:     for (j=0; j<3 && i>0; j++) {
551:       /* h <- A[0:i-1]'*m */
552:       SlepcMatDenseMult(h,i,0.0,1.0,A,ld,n,i,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
553:       /* h <- diag(ns)*h */
554:       for (k=0; k<i; k++) h[k] *= ns_[k];
555:       /* A[i] <- A[i] - A[0:i-1]*h */
556:       SlepcMatDenseMult(&A[ld*i],ld,1.0,-1.0,A,ld,n,i,PETSC_FALSE,h,i,i,1,PETSC_FALSE);
557:       /* m <- diag(s)*A[i] */
558:       for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
559:       /* nr_o <- mynorm(A[i]'*m) */
560:       SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
561:       nr = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
562:       if (PetscAbs(nr) < PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Linear dependency detected");
563:       if (PetscAbs(nr) > 0.7*PetscAbs(nr_o)) break;
564:       nr_o = nr;
565:     }
566:     ns_[i] = PetscSign(nr);
567:     /* A[i] <- A[i]/|nr| */
568:     nr_abs = PetscAbs(nr);
569:     PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nr_abs,&done,&rA_,&one,A+i*ld,&ld_,&info));
570:     SlepcCheckLapackInfo("lascl",info);
571:   }
572:   PetscLogEventEnd(DS_Other,ds,0,0,0);
573:   DSRestoreArray(ds,mat,&A_);
574:   PetscObjectStateIncrease((PetscObject)ds);
575:   if (lindcols) *lindcols = cols;
576:   return(0);
577: }