~ubuntu-branches/debian/sid/arpack/sid : revision 14

1

c\BeginDoc

2

c

3

c\Name: pdseupd

3

c\Name: pdseupd

4

c

5

c Message Passing Layer: MPI

5

c Message Passing Layer: MPI

6

c

7

c\Description:

8

c

41

c with a single call.

42

c

43

c\Usage:

44

c call pdseupd

44

c call pdseupd

45

c ( COMM, RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, BMAT, N, WHICH, NEV, TOL,

46

c RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD, WORKL, LWORKL, INFO )

47

c

48

c\Arguments

49

c COMM MPI Communicator for the processor grid. (INPUT)

49

c COMM MPI Communicator for the processor grid. (INPUT)

50

c

51

c RVEC LOGICAL (INPUT)

52

c Specifies whether Ritz vectors corresponding to the Ritz value

69

c Ritz value D(j), SELECT(j) must be set to .TRUE..

70

c If HOWMNY = 'A' , SELECT is used as workspace.

71

c

72

c D Double precision array of dimension NEV. (OUTPUT)

72

c D Double precision array of dimension NEV. (OUTPUT)

73

c On exit, D contains the Ritz value approximations to the

74

c eigenvalues of A*z = lambda*B*z. The values are returned

75

c in ascending order. If IPARAM(7) = 3,4,5 then D represents

76

c the Ritz values of OP computed by pdsaupd transformed to

76

c the Ritz values of OP computed by pdsaupd transformed to

77

c those of the original eigensystem A*z = lambda*B*z. If

78

c IPARAM(7) = 1,2 then the Ritz values of OP are the same

79

c as the those of A*z = lambda*B*z.

80

c

81

c Z Double precision N by NEV array if HOWMNY = 'A'. (OUTPUT)

81

c Z Double precision N by NEV array if HOWMNY = 'A'. (OUTPUT)

82

c On exit, Z contains the B-orthonormal Ritz vectors of the

83

c eigensystem A*z = lambda*B*z corresponding to the Ritz

84

c value approximations.

90

c The leading dimension of the array Z. If Ritz vectors are

91

c desired, then LDZ .ge. max( 1, N ). In any case, LDZ .ge. 1.

92

c

93

c SIGMA Double precision (INPUT)

93

c SIGMA Double precision (INPUT)

94

c If IPARAM(7) = 3,4,5 represents the shift. Not referenced if

95

c IPARAM(7) = 1 or 2.

96

c

97

c

98

c **** The remaining arguments MUST be the same as for the ****

99

c **** call to PDNAUPD that was just completed. ****

99

c **** call to PDNAUPD that was just completed. ****

100

c

101

c NOTE: The remaining arguments

102

c

109

c

110

c Two of these parameters (WORKL, INFO) are also output parameters:

111

c

112

c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)

112

c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)

113

c WORKL(1:4*ncv) contains information obtained in

114

c PSSAUPD. They are not changed by PSSEUPD.

115

c WORKL(4*ncv+1:ncv*ncv+8*ncv) holds the

136

c = -6: BMAT must be one of 'I' or 'G'.

137

c = -7: Length of private work WORKL array is not sufficient.

138

c = -8: Error return from trid. eigenvalue calculation;

139

c Information error from LAPACK routine dsteqr .

139

c Information error from LAPACK routine dsteqr.

140

c = -9: Starting vector is zero.

141

c = -10: IPARAM(7) must be 1,2,3,4,5.

142

c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.

145

c accuracy.

146

c = -15: HOWMNY must be one of 'A' or 'S' if RVEC = .true.

147

c = -16: HOWMNY = 'S' not yet implemented

148

c = -17: DSEUPD got a different count of the number of converged

149

c Ritz values than DSAUPD got. This indicates the user

150

c probably made an error in passing data from DSAUPD to

151

c DSEUPD or that the data was modified before entering

152

c DSEUPD .

148

c = -17: DSEUPD got a different count of the number of converged

149

c Ritz values than DSAUPD got. This indicates the user

150

c probably made an error in passing data from DSAUPD to

151

c DSEUPD or that the data was modified before entering

152

c DSEUPD.

153

c

154

c\BeginLib

155

c

182

c stage for the user who wants to incorporate it.

183

c

184

c\Routines called:

185

c dsesrt ARPACK routine that sorts an array X, and applies the

185

c dsesrt ARPACK routine that sorts an array X, and applies the

186

c corresponding permutation to a matrix A.

187

c dsortr dsortr ARPACK sorting routine.

188

c pdnorm2 Parallel ARPACK routine that computes the 2-norm of a vector.

187

c dsortr dsortr ARPACK sorting routine.

188

c pdnorm2 Parallel ARPACK routine that computes the 2-norm of a vector.

189

c pivout Parallel ARPACK utility routine that prints integers.

190

c pdvout Parallel ARPACK utility routine that prints vectors.

191

c dgeqr2 LAPACK routine that computes the QR factorization of

190

c pdvout Parallel ARPACK utility routine that prints vectors.

191

c dgeqr2 LAPACK routine that computes the QR factorization of

192

c a matrix.

193

c dlacpy LAPACK matrix copy routine.

194

c pdlamch ScaLAPACK routine that determines machine constants.

195

c dorm2r LAPACK routine that applies an orthogonal matrix in

193

c dlacpy LAPACK matrix copy routine.

194

c pdlamch ScaLAPACK routine that determines machine constants.

195

c dorm2r LAPACK routine that applies an orthogonal matrix in

196

c factored form.

197

c dsteqr LAPACK routine that computes eigenvalues and eigenvectors

197

c dsteqr LAPACK routine that computes eigenvalues and eigenvectors

198

c of a tridiagonal matrix.

199

c dger Level 2 BLAS rank one update to a matrix.

200

c dcopy Level 1 BLAS that copies one vector to another .

201

c dscal Level 1 BLAS that scales a vector.

202

c dswap Level 1 BLAS that swaps the contents of two vectors.

199

c dger Level 2 BLAS rank one update to a matrix.

200

c dcopy Level 1 BLAS that copies one vector to another .

201

c dscal Level 1 BLAS that scales a vector.

202

c dswap Level 1 BLAS that swaps the contents of two vectors.

203

c\Authors

204

c Danny Sorensen Phuong Vu

205

c Richard Lehoucq CRPC / Rice University

216

c Starting Point: Serial Code FILE: seupd.F SID: 2.4

217

c

218

c\SCCS Information:

219

c FILE: seupd.F SID: 1.10 DATE OF SID: 04/10/01

219

c FILE: seupd.F SID: 1.11 DATE OF SID: 10/25/03

220

c

221

c\EndLib

222

c

223

c-----------------------------------------------------------------------

224

subroutine pdseupd

224

subroutine pdseupd

225

& (comm , rvec , howmny, select, d ,

226

& z , ldz , sigma , bmat , n ,

227

& which , nev , tol , resid , ncv ,

229

& workl , lworkl, info )

230

c

231

c %--------------------%

232

c | MPI Communicator |

232

c | MPI Communicator |

233

c %--------------------%

234

c

235

integer comm

248

character bmat, howmny, which*2

249

logical rvec

250

integer info, ldz, ldv, lworkl, n, ncv, nev

251

Double precision

251

Double precision

252

& sigma, tol

253

c

254

c %-----------------%

257

c

258

integer iparam(7), ipntr(11)

259

logical select(ncv)

260

Double precision

260

Double precision

261

& d(nev), resid(n), v(ldv,ncv), z(ldz, nev),

262

& workd(2*n), workl(lworkl)

263

c

265

c | Parameters |

266

c %------------%

267

c

268

Double precision

268

Double precision

269

& one, zero

270

parameter (one = 1.0 , zero = 0.0 )

270

parameter (one = 1.0, zero = 0.0)

271

c

272

c %---------------%

273

c | Local Scalars |

279

& ldq , mode , msglvl, nconv , next ,

280

& ritz , irz , ibd , np , ishift,

281

& leftptr, rghtptr, numcnv, jj

282

Double precision

282

Double precision

283

& bnorm2, rnorm, temp, temp1, eps23

284

logical reord

285

c

287

c | External Subroutines |

288

c %----------------------%

289

c

290

external dcopy , dger , dgeqr2 , dlacpy , dorm2r , dscal ,

291

& dsesrt , dsteqr , dswap , pdvout , pivout, dsortr

290

external dcopy , dger , dgeqr2, dlacpy, dorm2r, dscal,

291

& dsesrt, dsteqr, dswap , pdvout, pivout, dsortr

292

c

293

c %--------------------%

294

c | External Functions |

295

c %--------------------%

296

c

297

Double precision

298

& pdnorm2 , pdlamch

299

external pdnorm2 , pdlamch

297

Double precision

298

& pdnorm2, pdlamch

299

external pdnorm2, pdlamch

300

c

301

c %---------------------%

302

c | Intrinsic Functions |

373

c | workl(1:2*ncv) := generated tridiagonal matrix H |

374

c | The subdiagonal is stored in workl(2:ncv). |

375

c | The dead spot is workl(1) but upon exiting |

376

c | pdsaupd stores the B-norm of the last residual |

376

c | pdsaupd stores the B-norm of the last residual |

377

c | vector in workl(1). We use this !!! |

378

c | workl(2*ncv+1:2*ncv+ncv) := ritz values |

379

c | The wanted values are in the first NCONV spots. |

380

c | workl(3*ncv+1:3*ncv+ncv) := computed Ritz estimates |

381

c | The wanted values are in the first NCONV spots. |

382

c | NOTE: workl(1:4*ncv) is set by pdsaupd and is not |

383

c | modified by pdseupd . |

382

c | NOTE: workl(1:4*ncv) is set by pdsaupd and is not |

383

c | modified by pdseupd. |

384

c %-------------------------------------------------------%

385

c

386

c %-------------------------------------------------------%

387

c | The following is used and set by pdseupd . |

387

c | The following is used and set by pdseupd. |

388

c | workl(4*ncv+1:4*ncv+ncv) := used as workspace during |

389

c | computation of the eigenvectors of H. Stores |

390

c | the diagonal of H. Upon EXIT contains the NCV |

400

c | workl(3*ncv+1:4*ncv). |

401

c | workl(6*ncv+1:6*ncv+ncv*ncv) := orthogonal Q that is |

402

c | the eigenvector matrix for H as returned by |

403

c | dsteqr . Not referenced if RVEC = .False. |

403

c | dsteqr. Not referenced if RVEC = .False. |

404

c | Ordering follows that of workl(4*ncv+1:5*ncv) |

405

c | workl(6*ncv+ncv*ncv+1:6*ncv+ncv*ncv+2*ncv) := |

406

c | Workspace. Needed by dsteqr and by pdseupd . |

406

c | Workspace. Needed by dsteqr and by pdseupd. |

407

c | GRAND total of NCV*(NCV+8) locations. |

408

c %-------------------------------------------------------%

409

c

439

c | Set machine dependent constant. |

440

c %---------------------------------%

441

c

442

eps23 = pdlamch (comm, 'Epsilon-Machine')

443

eps23 = eps23**(2.0 / 3.0 )

442

eps23 = pdlamch(comm, 'Epsilon-Machine')

443

eps23 = eps23**(2.0 / 3.0)

444

c

445

c %---------------------------------------%

446

c | RNORM is B-norm of the RESID(1:N). |

447

c | BNORM2 is the 2 norm of B*RESID(1:N). |

448

c | Upon exit of pdsaupd WORKD(1:N) has |

448

c | Upon exit of pdsaupd WORKD(1:N) has |

449

c | B*RESID(1:N). |

450

c %---------------------------------------%

451

c

453

if (bmat .eq. 'I') then

454

bnorm2 = rnorm

455

else if (bmat .eq. 'G') then

456

bnorm2 = pdnorm2 (comm, n, workd, 1)

456

bnorm2 = pdnorm2(comm, n, workd, 1)

457

end if

458

c

459

if (msglvl .gt. 2) then

460

call pdvout (comm, logfil, ncv, workl(irz), ndigit,

460

call pdvout(comm, logfil, ncv, workl(irz), ndigit,

461

& '_seupd: Ritz values passed in from _SAUPD.')

462

call pdvout (comm, logfil, ncv, workl(ibd), ndigit,

462

call pdvout(comm, logfil, ncv, workl(ibd), ndigit,

463

& '_seupd: Ritz estimates passed in from _SAUPD.')

464

end if

465

if (rvec) then

487

c

488

np = ncv - nev

489

ishift = 0

490

call pdsgets (comm , ishift, which ,

490

call pdsgets(comm , ishift, which ,

491

& nev , np , workl(irz),

492

& workl(bounds), workl , workl(np+1))

492

& workl(bounds), workl)

493

c

494

if (msglvl .gt. 2) then

495

call pdvout (comm, logfil, ncv, workl(irz), ndigit,

495

call pdvout(comm, logfil, ncv, workl(irz), ndigit,

496

& '_seupd: Ritz values after calling _SGETS.')

497

call pdvout (comm, logfil, ncv, workl(bounds), ndigit,

497

call pdvout(comm, logfil, ncv, workl(bounds), ndigit,

498

& '_seupd: Ritz value indices after calling _SGETS.')

499

end if

500

c

540

c | Initialize the eigenvector matrix Q to the identity. |

541

c %-----------------------------------------------------------%

542

c

543

call dcopy (ncv-1, workl(ih+1) , 1, workl(ihb), 1)

544

call dcopy (ncv , workl(ih+ldh), 1, workl(ihd), 1)

543

call dcopy (ncv-1, workl(ih+1) , 1, workl(ihb), 1)

544

call dcopy (ncv , workl(ih+ldh), 1, workl(ihd), 1)

545

c

546

call dsteqr ('Identity', ncv , workl(ihd),

546

call dsteqr('Identity', ncv , workl(ihd),

547

& workl(ihb), workl(iq), ldq ,

548

& workl(iw) , ierr)

549

c

553

end if

554

c

555

if (msglvl .gt. 1) then

556

call dcopy (ncv, workl(iq+ncv-1), ldq, workl(iw), 1)

557

call pdvout (comm, logfil, ncv, workl(ihd), ndigit,

556

call dcopy (ncv, workl(iq+ncv-1), ldq, workl(iw), 1)

557

call pdvout (comm, logfil, ncv, workl(ihd), ndigit,

558

& '_seupd: NCV Ritz values of the final H matrix')

559

call pdvout (comm, logfil, ncv, workl(iw), ndigit,

559

call pdvout (comm, logfil, ncv, workl(iw), ndigit,

560

& '_seupd: last row of the eigenvector matrix for H')

561

end if

562

c

607

temp = workl(ihd+leftptr-1)

608

workl(ihd+leftptr-1) = workl(ihd+rghtptr-1)

609

workl(ihd+rghtptr-1) = temp

610

call dcopy (ncv, workl(iq+ncv*(leftptr-1)), 1,

610

call dcopy(ncv, workl(iq+ncv*(leftptr-1)), 1,

611

& workl(iw), 1)

612

call dcopy (ncv, workl(iq+ncv*(rghtptr-1)), 1,

612

call dcopy(ncv, workl(iq+ncv*(rghtptr-1)), 1,

613

& workl(iq+ncv*(leftptr-1)), 1)

614

call dcopy (ncv, workl(iw), 1,

614

call dcopy(ncv, workl(iw), 1,

615

& workl(iq+ncv*(rghtptr-1)), 1)

616

leftptr = leftptr + 1

617

rghtptr = rghtptr - 1

623

30 end if

624

c

625

if (msglvl .gt. 2) then

626

call pdvout (comm, logfil, ncv, workl(ihd), ndigit,

626

call pdvout (comm, logfil, ncv, workl(ihd), ndigit,

627

& '_seupd: The eigenvalues of H--reordered')

628

end if

629

c

631

c | Load the converged Ritz values into D. |

632

c %----------------------------------------%

633

c

634

call dcopy (nconv, workl(ihd), 1, d, 1)

634

call dcopy(nconv, workl(ihd), 1, d, 1)

635

c

636

else

637

c

639

c | Ritz vectors not required. Load Ritz values into D. |

640

c %-----------------------------------------------------%

641

c

642

call dcopy (nconv, workl(ritz), 1, d, 1)

643

call dcopy (ncv, workl(ritz), 1, workl(ihd), 1)

642

call dcopy(nconv, workl(ritz), 1, d, 1)

643

call dcopy(ncv, workl(ritz), 1, workl(ihd), 1)

644

c

645

end if

646

c

658

c %---------------------------------------------------------%

659

c

660

if (rvec) then

661

call dsesrt ('LA', rvec , nconv, d, ncv, workl(iq), ldq)

661

call dsesrt('LA', rvec , nconv, d, ncv, workl(iq), ldq)

662

else

663

call dcopy (ncv, workl(bounds), 1, workl(ihb), 1)

663

call dcopy(ncv, workl(bounds), 1, workl(ihb), 1)

664

end if

665

c

666

else

674

c | lambda = sigma * theta / ( theta - 1 ) |

675

c | For TYPE = 'CAYLEY' the transformation is |

676

c | lambda = sigma * (theta + 1) / (theta - 1 ) |

677

c | where the theta are the Ritz values returned by pdsaupd . |

677

c | where the theta are the Ritz values returned by pdsaupd. |

678

c | NOTES: |

679

c | *The Ritz vectors are not affected by the transformation. |

680

c | They are only reordered. |

681

c %-------------------------------------------------------------%

682

c

683

call dcopy (ncv, workl(ihd), 1, workl(iw), 1)

683

call dcopy (ncv, workl(ihd), 1, workl(iw), 1)

684

if (type .eq. 'SHIFTI') then

685

do 40 k=1, ncv

686

workl(ihd+k-1) = one / workl(ihd+k-1) + sigma

712

c | Ritz vector purification. |

713

c %-------------------------------------------------------------%

714

c

715

call dcopy (nconv, workl(ihd), 1, d, 1)

716

call dsortr ('LA', .true., nconv, workl(ihd), workl(iw))

715

call dcopy (nconv, workl(ihd), 1, d, 1)

716

call dsortr('LA', .true., nconv, workl(ihd), workl(iw))

717

if (rvec) then

718

call dsesrt ('LA', rvec , nconv, d, ncv, workl(iq), ldq)

718

call dsesrt('LA', rvec , nconv, d, ncv, workl(iq), ldq)

719

else

720

call dcopy (ncv, workl(bounds), 1, workl(ihb), 1)

721

call dscal (ncv, bnorm2/rnorm, workl(ihb), 1)

722

call dsortr ('LA', .true., nconv, d, workl(ihb))

720

call dcopy(ncv, workl(bounds), 1, workl(ihb), 1)

721

call dscal(ncv, bnorm2/rnorm, workl(ihb), 1)

722

call dsortr('LA', .true., nconv, d, workl(ihb))

723

end if

724

c

725

end if

738

c | columns of workl(iq,ldq). |

739

c %----------------------------------------------------------%

740

c

741

call dgeqr2 (ncv, nconv , workl(iq) ,

741

call dgeqr2(ncv, nconv , workl(iq) ,

742

& ldq, workl(iw+ncv), workl(ihb),

743

& ierr)

744

c

750

c | the Ritz values in workl(ihd). |

751

c %--------------------------------------------------------%

752

c

753

call dorm2r ('Right' , 'Notranspose', n ,

753

call dorm2r('Right' , 'Notranspose', n ,

754

& ncv , nconv , workl(iq),

755

& ldq , workl(iw+ncv), v ,

756

& ldv , workd(n+1) , ierr )

757

call dlacpy ('All', n, nconv, v, ldv, z, ldz)

757

call dlacpy('All', n, nconv, v, ldv, z, ldz)

758

c

759

c %-----------------------------------------------------%

760

c | In order to compute the Ritz estimates for the Ritz |

766

workl(ihb+j-1) = zero

767

65 continue

768

workl(ihb+ncv-1) = one

769

call dorm2r ('Left', 'Transpose' , ncv ,

769

call dorm2r('Left', 'Transpose' , ncv ,

770

& 1 , nconv , workl(iq) ,

771

& ldq , workl(iw+ncv), workl(ihb),

772

& ncv , temp , ierr )

790

c | If RVEC = .true. then compute Ritz estimates |

791

c | of the theta. |

792

c | If RVEC = .false. then copy Ritz estimates |

793

c | as computed by pdsaupd . |

793

c | as computed by pdsaupd. |

794

c | * Determine Ritz estimates of the lambda. |

795

c %-------------------------------------------------%

796

c

797

call dscal (ncv, bnorm2, workl(ihb), 1)

797

call dscal (ncv, bnorm2, workl(ihb), 1)

798

if (type .eq. 'SHIFTI') then

799

c

800

do 80 k=1, ncv

821

end if

822

c

823

if (type .ne. 'REGULR' .and. msglvl .gt. 1) then

824

call pdvout (comm, logfil, nconv, d, ndigit,

824

call pdvout (comm, logfil, nconv, d, ndigit,

825

& '_seupd: Untransformed converged Ritz values')

826

call pdvout (comm, logfil, nconv, workl(ihb), ndigit,

826

call pdvout (comm, logfil, nconv, workl(ihb), ndigit,

827

& '_seupd: Ritz estimates of the untransformed Ritz values')

828

else if (msglvl .gt. 1) then

829

call pdvout (comm, logfil, nconv, d, ndigit,

829

call pdvout (comm, logfil, nconv, d, ndigit,

830

& '_seupd: Converged Ritz values')

831

call pdvout (comm, logfil, nconv, workl(ihb), ndigit,

831

call pdvout (comm, logfil, nconv, workl(ihb), ndigit,

832

& '_seupd: Associated Ritz estimates')

833

end if

834

c

855

end if

856

c

857

if (type .ne. 'REGULR')

858

& call dger (n, nconv, one, resid, 1, workl(iw), 1, z, ldz)

858

& call dger(n, nconv, one, resid, 1, workl(iw), 1, z, ldz)

859

c

860

9000 continue

861

c

862

return

863

c

864

c %----------------%

865

c | End of pdseupd |

865

c | End of pdseupd |

866

c %----------------%

867

c

868

end