Compare commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Sylvestre Ledru
|
9e7ceb7480 | ||
|
Sylvestre Ledru
|
e7ad0850a1 | ||
|
Franck HOUSSEN
|
ee9a93a853 | ||
|
Franck HOUSSEN
|
2d969e902d | ||
|
Franck HOUSSEN
|
5872987e67 | ||
|
Franck HOUSSEN
|
bbb87ba0ff | ||
|
Franck HOUSSEN
|
78d781286a | ||
|
Franck HOUSSEN
|
3829c449c0 | ||
|
Franck HOUSSEN
|
b7ce7debe5 | ||
|
caliarim
|
31854cadaf |
@@ -40,6 +40,8 @@ TESTS/bug_58_double
|
||||
TESTS/bug_79_double_complex
|
||||
TESTS/icb_arpack_c
|
||||
TESTS/icb_arpack_cpp
|
||||
TESTS/bug_142
|
||||
TESTS/bug_142_gen
|
||||
EXAMPLES/BAND/[sd][sn]bdr[123456]
|
||||
EXAMPLES/BAND/[cz]nbdr[1234]
|
||||
EXAMPLES/COMPLEX/[cz]ndrv[1234]
|
||||
|
||||
@@ -1,3 +1,13 @@
|
||||
arpack-ng - 3.6.3
|
||||
|
||||
[ Franck Houssen ]
|
||||
* Add Fortran common initialization (block data).
|
||||
|
||||
[ Marco Caliari ]
|
||||
* Give up forcing the initial residual to be in the range of the operator OP after a restart (Closes: #142).
|
||||
|
||||
-- Sylvestre Ledru <sylvestre@debian.org> Wed, 19 Sep 2018 09:59:59 +0200
|
||||
|
||||
arpack-ng - 3.6.2
|
||||
|
||||
* Remove all trailing whitespaces
|
||||
@@ -6,7 +16,7 @@ arpack-ng - 3.6.2
|
||||
* Install: move headers into a dedicated directory (local/include/arpack).
|
||||
(Closes #126)
|
||||
* Add configuration summary.
|
||||
* Improve the flag detection. Hopefully fix the ppc64el and other archs
|
||||
* Improve the flag detection. Hopefully fix the ppc64el and other archs
|
||||
issues in Debian
|
||||
|
||||
-- Sylvestre Ledru <sylvestre@debian.org> Sat, 23 Jun 2018 14:56:54 +0200
|
||||
|
||||
+12
-4
@@ -174,13 +174,13 @@ endif()
|
||||
# file(GLOB arpacklapack_STAT_SRCS lapack/*.f)
|
||||
#endif()
|
||||
#file(GLOB arpackutil_STAT_SRCS util/*.f)
|
||||
file(GLOB arpacksrc_STAT_SRCS ${arpack_SOURCE_DIR}/SRC/*.f)
|
||||
file(GLOB arpacksrc_STAT_SRCS dbgini.f staini.f ${arpack_SOURCE_DIR}/SRC/*.f)
|
||||
|
||||
set(arpacksrc_ICB "")
|
||||
set(parpacksrc_ICB "")
|
||||
if(ICB)
|
||||
file(GLOB arpacksrc_ICB SRC/icba*.f90 debug_init.f90 debug_icb.f90 stat_icb.f90)
|
||||
file(GLOB parpacksrc_ICB PARPACK/SRC/MPI/icbp*.f90 debug_init.f90 debug_icb.f90 stat_icb.f90)
|
||||
file(GLOB arpacksrc_ICB SRC/icba*.f90 debug_icb.f90 stat_icb.f90)
|
||||
file(GLOB parpacksrc_ICB PARPACK/SRC/MPI/icbp*.f90 debug_icb.f90 stat_icb.f90)
|
||||
endif()
|
||||
|
||||
set(arpackutil_STAT_SRCS
|
||||
@@ -202,7 +202,7 @@ set(arpackutil_STAT_SRCS
|
||||
|
||||
if (MPI)
|
||||
file(GLOB parpackutil_STAT_SRCS PARPACK/UTIL/MPI/*.f)
|
||||
file(GLOB parpacksrc_STAT_SRCS PARPACK/SRC/MPI/*.f)
|
||||
file(GLOB parpacksrc_STAT_SRCS dbgini.f staini.f PARPACK/SRC/MPI/*.f)
|
||||
endif()
|
||||
|
||||
# use -DBUILD_SHARED_LIBS=ON|OFF to control static/shared
|
||||
@@ -483,6 +483,14 @@ add_executable(bug_79_double_complex TESTS/bug_79_double_complex.f)
|
||||
target_link_libraries(bug_79_double_complex arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
|
||||
add_test(bug_79_double_complex_tst Tests/bug_79_double_complex)
|
||||
|
||||
add_executable(bug_142 TESTS/bug_142.f)
|
||||
target_link_libraries(bug_142 arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
|
||||
add_test(bug_142_tst Tests/bug_142)
|
||||
|
||||
add_executable(bug_142_gen TESTS/bug_142_gen.f)
|
||||
target_link_libraries(bug_142_gen arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
|
||||
add_test(bug_142_gen_tst Tests/bug_142_gen)
|
||||
|
||||
if(MPI)
|
||||
add_executable(issue46 PARPACK/TESTS/MPI/issue46.f)
|
||||
target_link_libraries(issue46 parpack arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
|
||||
|
||||
@@ -9,8 +9,12 @@ CNDRV = pcndrv1
|
||||
ZNDRV = pzndrv1
|
||||
|
||||
check_PROGRAMS = $(SNDRV) $(DNDRV) $(SSDRV) $(DSDRV) $(CNDRV) $(ZNDRV)
|
||||
# Disable tests for now. Needs to be run using mpirun
|
||||
#TESTS = $(check_PROGRAMS)
|
||||
|
||||
# Run MPI tests with "mpirun -n 2"
|
||||
LOG_COMPILER = mpirun
|
||||
LOG_FLAGS = -n 2
|
||||
|
||||
TESTS = $(check_PROGRAMS)
|
||||
|
||||
EXTRA_DIST = debug.h stat.h
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
AUTOMAKE_OPTIONS = subdir-objects # Needed as debug_init/icb.f90 are not in current directory.
|
||||
AUTOMAKE_OPTIONS = subdir-objects # Needed as debug/stat* are not in current directory.
|
||||
|
||||
F77 = $(MPIF77)
|
||||
FFLAGS_SAV = @FFLAGS@
|
||||
@@ -34,7 +34,8 @@ libparpack@LIBSUFFIX@_noopt_la_FFLAGS = -O0
|
||||
|
||||
lib_LTLIBRARIES = libparpack@LIBSUFFIX@.la
|
||||
libparpack@LIBSUFFIX@_la_SOURCES = $(PSRC) $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
|
||||
libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_init.f90
|
||||
libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/dbgini.f
|
||||
libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/staini.f
|
||||
if ICB
|
||||
libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_icb.f90
|
||||
libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/stat_icb.f90
|
||||
|
||||
@@ -120,7 +120,7 @@ int zn() {
|
||||
int select[ncv];
|
||||
double _Complex z[(N+1)*(nev+1)];
|
||||
BLASINT ldz = N+1;
|
||||
double sigma=0;
|
||||
double _Complex sigma=0. + I*0.;
|
||||
int k;
|
||||
for (k=0; k < 3*N; ++k )
|
||||
workd[k] = 0. + I * 0.;
|
||||
|
||||
@@ -113,7 +113,7 @@ void complex_symmetric_runner() {
|
||||
|
||||
float tol = 0.0f;
|
||||
bool rvec = true;
|
||||
float sigma = 0.0f;
|
||||
std::complex<float> sigma(0.0f, 0.0f);
|
||||
|
||||
std::vector<std::complex<float>> resid(N);
|
||||
std::vector<std::complex<float>> V(ncv * N);
|
||||
|
||||
@@ -255,7 +255,7 @@ c %-------------------------------------------%
|
||||
c
|
||||
rvec = .true.
|
||||
c
|
||||
call pdseupd ( comm, rvec, 'All', select,
|
||||
call pdseupd ( comm, rvec, 'A', select,
|
||||
& d, v, ldv, sigma,
|
||||
& bmat, nloc, which, nev, tol, resid, ncv, v, ldv,
|
||||
& iparam, ipntr, workd, workl, lworkl, ierr )
|
||||
|
||||
+3
-2
@@ -1,4 +1,4 @@
|
||||
AUTOMAKE_OPTIONS = subdir-objects # Needed as debug_init/icb.f90 are not in current directory.
|
||||
AUTOMAKE_OPTIONS = subdir-objects # Needed as debug/stat* are not in current directory.
|
||||
|
||||
SSRC = snaitr.f snapps.f snaup2.f snaupd.f snconv.f sneigh.f sneupd.f sngets.f sstatn.f \
|
||||
ssaitr.f ssapps.f ssaup2.f ssaupd.f ssconv.f sseigt.f sseupd.f ssgets.f sstats.f \
|
||||
@@ -25,7 +25,8 @@ EXTRA_DIST = debug.h stat.h version.h
|
||||
|
||||
lib_LTLIBRARIES = libarpack@LIBSUFFIX@.la
|
||||
libarpack@LIBSUFFIX@_la_SOURCES = $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
|
||||
libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_init.f90
|
||||
libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/dbgini.f
|
||||
libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/staini.f
|
||||
if ICB
|
||||
libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_icb.f90
|
||||
libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/stat_icb.f90
|
||||
|
||||
+11
-7
@@ -239,12 +239,16 @@ c | the generalized problem when B is possibly (singular). |
|
||||
c %----------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call ccopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
if (itry .eq. 1) then
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call ccopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
else if (itry .gt. 1 .and. bmat .eq. 'G') then
|
||||
call ccopy (n, resid, 1, workd(n + 1), 1)
|
||||
end if
|
||||
end if
|
||||
c
|
||||
c %----------------------------------------%
|
||||
@@ -269,7 +273,7 @@ c %------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
first = .TRUE.
|
||||
call ccopy (n, workd(n+1), 1, resid, 1)
|
||||
if (itry .eq. 1) call ccopy (n, workd(n + 1), 1, resid, 1)
|
||||
if (bmat .eq. 'G') then
|
||||
nbx = nbx + 1
|
||||
ipntr(1) = n + 1
|
||||
|
||||
+11
-7
@@ -242,12 +242,16 @@ c | the generalized problem when B is possibly (singular). |
|
||||
c %----------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call dcopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
if (itry .eq. 1) then
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call dcopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
else if (itry .gt. 1 .and. bmat .eq. 'G') then
|
||||
call dcopy (n, resid, 1, workd(n + 1), 1)
|
||||
end if
|
||||
end if
|
||||
c
|
||||
c %-----------------------------------------%
|
||||
@@ -274,7 +278,7 @@ c %------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
first = .TRUE.
|
||||
call dcopy (n, workd(n+1), 1, resid, 1)
|
||||
if (itry .eq. 1) call dcopy (n, workd(n + 1), 1, resid, 1)
|
||||
if (bmat .eq. 'G') then
|
||||
nbx = nbx + 1
|
||||
ipntr(1) = n + 1
|
||||
|
||||
+11
-7
@@ -242,12 +242,16 @@ c | the generalized problem when B is possibly (singular). |
|
||||
c %----------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call scopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
if (itry .eq. 1) then
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call scopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
else if (itry .gt. 1 .and. bmat .eq. 'G') then
|
||||
call scopy (n, resid, 1, workd(n + 1), 1)
|
||||
end if
|
||||
end if
|
||||
c
|
||||
c %-----------------------------------------%
|
||||
@@ -274,7 +278,7 @@ c %------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
first = .TRUE.
|
||||
call scopy (n, workd(n+1), 1, resid, 1)
|
||||
if (itry .eq. 1) call scopy (n, workd(n + 1), 1, resid, 1)
|
||||
if (bmat .eq. 'G') then
|
||||
nbx = nbx + 1
|
||||
ipntr(1) = n + 1
|
||||
|
||||
+11
-7
@@ -239,12 +239,16 @@ c | the generalized problem when B is possibly (singular). |
|
||||
c %----------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call zcopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
if (itry .eq. 1) then
|
||||
nopx = nopx + 1
|
||||
ipntr(1) = 1
|
||||
ipntr(2) = n + 1
|
||||
call zcopy (n, resid, 1, workd, 1)
|
||||
ido = -1
|
||||
go to 9000
|
||||
else if (itry .gt. 1 .and. bmat .eq. 'G') then
|
||||
call zcopy (n, resid, 1, workd(n + 1), 1)
|
||||
end if
|
||||
end if
|
||||
c
|
||||
c %----------------------------------------%
|
||||
@@ -269,7 +273,7 @@ c %------------------------------------------------------%
|
||||
c
|
||||
call arscnd (t2)
|
||||
first = .TRUE.
|
||||
call zcopy (n, workd(n+1), 1, resid, 1)
|
||||
if (itry .eq. 1) call zcopy (n, workd(n + 1), 1, resid, 1)
|
||||
if (bmat .eq. 'G') then
|
||||
nbx = nbx + 1
|
||||
ipntr(1) = n + 1
|
||||
|
||||
+5
-1
@@ -1,6 +1,6 @@
|
||||
LDADD=$(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
|
||||
|
||||
check_PROGRAMS = dnsimp bug_1323 bug_58_double bug_79_double_complex
|
||||
check_PROGRAMS = dnsimp bug_1323 bug_58_double bug_79_double_complex bug_142 bug_142_gen
|
||||
if ICB
|
||||
check_PROGRAMS += icb_arpack_c
|
||||
check_PROGRAMS += icb_arpack_cpp
|
||||
@@ -21,6 +21,10 @@ bug_58_double_SOURCES = bug_58_double.f
|
||||
|
||||
bug_79_double_complex_SOURCES = bug_79_double_complex.f
|
||||
|
||||
bug_142_SOURCES = bug_142.f
|
||||
|
||||
bug_142_gen_SOURCES = bug_142_gen.f
|
||||
|
||||
if ICB
|
||||
icb_arpack_c_SOURCES = icb_arpack_c.c
|
||||
icb_arpack_c_LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la
|
||||
|
||||
+414
@@ -0,0 +1,414 @@
|
||||
program bug_142
|
||||
c
|
||||
c Avoid taking the initial vector in the range of OP after a restart
|
||||
c (standard case)
|
||||
c
|
||||
c Example program to illustrate the idea of reverse communication
|
||||
c for a standard nonsymmetric eigenvalue problem.
|
||||
c
|
||||
c We implement example one of ex-nonsym.doc in DOCUMENTS directory
|
||||
c
|
||||
c\Example-1
|
||||
c ... Suppose we want to solve A*x = lambda*x in regular mode,
|
||||
c where A is obtained from the standard central difference
|
||||
c discretization of the convection-diffusion operator
|
||||
c (Laplacian u) + rho*(du / dx)
|
||||
c on the unit square [0,1]x[0,1] with zero Dirichlet boundary
|
||||
c condition.
|
||||
c
|
||||
c ... OP = A and B = I.
|
||||
c
|
||||
c ... Assume "call av (nx,x,y)" computes y = A*x.c
|
||||
c
|
||||
c ... Use mode 1 of DNAUPD.
|
||||
c
|
||||
c\BeginLib
|
||||
c
|
||||
c\Routines called:
|
||||
c dnaupd ARPACK reverse communication interface routine.
|
||||
c dneupd ARPACK routine that returns Ritz values and (optionally)
|
||||
c Ritz vectors.
|
||||
c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
|
||||
c daxpy Level 1 BLAS that computes y <- alpha*x+y.
|
||||
c dnrm2 Level 1 BLAS that computes the norm of a vector.
|
||||
c av Matrix vector multiplication routine that computes A*x.
|
||||
c tv Matrix vector multiplication routine that computes T*x,
|
||||
c where T is a tridiagonal matrix. It is used in routine
|
||||
c av.
|
||||
c
|
||||
c\Author
|
||||
c Richard Lehoucq
|
||||
c Danny Sorensen
|
||||
c Chao Yang
|
||||
c Dept. of Computational &
|
||||
c Applied Mathematics
|
||||
c Rice University
|
||||
c Houston, Texas
|
||||
c
|
||||
c\SCCS Information: @(#)
|
||||
c FILE: ndrv1.F SID: 2.5 DATE OF SID: 10/17/00 RELEASE: 2
|
||||
c
|
||||
c\Remarks
|
||||
c 1. None
|
||||
c
|
||||
c\EndLib
|
||||
c---------------------------------------------------------------------------
|
||||
c
|
||||
c %-----------------------------%
|
||||
c | Define maximum dimensions |
|
||||
c | for all arrays. |
|
||||
c | MAXN: Maximum dimension |
|
||||
c | of the A allowed. |
|
||||
c | MAXNEV: Maximum NEV allowed |
|
||||
c | MAXNCV: Maximum NCV allowed |
|
||||
c %-----------------------------%
|
||||
c
|
||||
integer maxn, maxnev, maxncv, ldv
|
||||
parameter (maxn=256, maxnev=12, maxncv=30, ldv=maxn)
|
||||
c
|
||||
c %--------------%
|
||||
c | Local Arrays |
|
||||
c %--------------%
|
||||
c
|
||||
integer iparam(11), ipntr(14)
|
||||
logical select(maxncv)
|
||||
Double precision
|
||||
& ax(maxn), d(maxncv,3), resid(maxn),
|
||||
& v(ldv,maxncv), workd(3*maxn),
|
||||
& workev(3*maxncv),
|
||||
& workl(3*maxncv*maxncv+6*maxncv), a(maxn, maxn)
|
||||
c
|
||||
c %---------------%
|
||||
c | Local Scalars |
|
||||
c %---------------%
|
||||
c
|
||||
character bmat*1, which*2
|
||||
integer ido, n, nx, nev, ncv, lworkl, info, j,
|
||||
& ierr, nconv, maxitr, ishfts, mode
|
||||
Double precision
|
||||
& tol, sigmar, sigmai
|
||||
logical first, rvec
|
||||
c
|
||||
c %------------%
|
||||
c | Parameters |
|
||||
c %------------%
|
||||
c
|
||||
Double precision
|
||||
& zero
|
||||
parameter (zero = 0.0D+0)
|
||||
c
|
||||
c %-----------------------------%
|
||||
c | BLAS & LAPACK routines used |
|
||||
c %-----------------------------%
|
||||
c
|
||||
Double precision
|
||||
& dlapy2, dnrm2
|
||||
external dlapy2, dnrm2, daxpy
|
||||
c
|
||||
c %--------------------%
|
||||
c | Intrinsic function |
|
||||
c %--------------------%
|
||||
c
|
||||
intrinsic abs
|
||||
c
|
||||
c %-----------------------%
|
||||
c | Executable Statements |
|
||||
c %-----------------------%
|
||||
c
|
||||
c %--------------------------------------------------%
|
||||
c | The number NX is the number of interior points |
|
||||
c | in the discretization of the 2-dimensional |
|
||||
c | convection-diffusion operator on the unit |
|
||||
c | square with zero Dirichlet boundary condition. |
|
||||
c | The number N(=NX*NX) is the dimension of the |
|
||||
c | matrix. A standard eigenvalue problem is |
|
||||
c | solved (BMAT = 'I'). NEV is the number of |
|
||||
c | eigenvalues to be approximated. The user can |
|
||||
c | modify NX, NEV, NCV, WHICH to solve problems of |
|
||||
c | different sizes, and to get different parts of |
|
||||
c | the spectrum. However, The following |
|
||||
c | conditions must be satisfied: |
|
||||
c | N <= MAXN |
|
||||
c | NEV <= MAXNEV |
|
||||
c | NEV + 2 <= NCV <= MAXNCV |
|
||||
c %--------------------------------------------------%
|
||||
c
|
||||
nx = 10
|
||||
n = 11
|
||||
nev = 1
|
||||
ncv = 11
|
||||
do i = 1,n
|
||||
do j = 1,n
|
||||
a(i,j) = 0.15d0/11
|
||||
end do
|
||||
end do
|
||||
do j = 2,n
|
||||
a(1,j) = a(1,j) + 0.85d0
|
||||
end do
|
||||
do i = 2,n
|
||||
a(i,1) = (1-a(1,1))/10
|
||||
end do
|
||||
|
||||
if ( n .gt. maxn ) then
|
||||
print *, ' ERROR with _NDRV1: N is greater than MAXN '
|
||||
go to 9000
|
||||
else if ( nev .gt. maxnev ) then
|
||||
print *, ' ERROR with _NDRV1: NEV is greater than MAXNEV '
|
||||
go to 9000
|
||||
else if ( ncv .gt. maxncv ) then
|
||||
print *, ' ERROR with _NDRV1: NCV is greater than MAXNCV '
|
||||
go to 9000
|
||||
end if
|
||||
bmat = 'I'
|
||||
which = 'LM'
|
||||
c
|
||||
c %-----------------------------------------------------%
|
||||
c | The work array WORKL is used in DNAUPD as |
|
||||
c | workspace. Its dimension LWORKL is set as |
|
||||
c | illustrated below. The parameter TOL determines |
|
||||
c | the stopping criterion. If TOL<=0, machine |
|
||||
c | precision is used. The variable IDO is used for |
|
||||
c | reverse communication, and is initially set to 0. |
|
||||
c | Setting INFO=0 indicates that a random vector is |
|
||||
c | generated in DNAUPD to start the Arnoldi iteration. |
|
||||
c %-----------------------------------------------------%
|
||||
c
|
||||
lworkl = 3*ncv**2+6*ncv
|
||||
tol = zero
|
||||
ido = 0
|
||||
info = 0
|
||||
c
|
||||
c %---------------------------------------------------%
|
||||
c | This program uses exact shifts with respect to |
|
||||
c | the current Hessenberg matrix (IPARAM(1) = 1). |
|
||||
c | IPARAM(3) specifies the maximum number of Arnoldi |
|
||||
c | iterations allowed. Mode 1 of DNAUPD is used |
|
||||
c | (IPARAM(7) = 1). All these options can be changed |
|
||||
c | by the user. For details see the documentation in |
|
||||
c | DNAUPD. |
|
||||
c %---------------------------------------------------%
|
||||
c
|
||||
ishfts = 1
|
||||
maxitr = 300
|
||||
mode = 1
|
||||
c
|
||||
iparam(1) = ishfts
|
||||
iparam(3) = maxitr
|
||||
iparam(7) = mode
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | M A I N L O O P (Reverse communication) |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
10 continue
|
||||
c
|
||||
c %---------------------------------------------%
|
||||
c | Repeatedly call the routine DNAUPD and take |
|
||||
c | actions indicated by parameter IDO until |
|
||||
c | either convergence is indicated or maxitr |
|
||||
c | has been exceeded. |
|
||||
c %---------------------------------------------%
|
||||
c
|
||||
call dnaupd ( ido, bmat, n, which, nev, tol, resid,
|
||||
& ncv, v, ldv, iparam, ipntr, workd, workl, lworkl,
|
||||
& info )
|
||||
c
|
||||
if (ido .eq. -1 .or. ido .eq. 1) then
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | Perform matrix vector multiplication |
|
||||
c | y <--- OP*x |
|
||||
c | The user should supply his/her own |
|
||||
c | matrix vector multiplication routine here |
|
||||
c | that takes workd(ipntr(1)) as the input |
|
||||
c | vector, and return the matrix vector |
|
||||
c | product to workd(ipntr(2)). |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& workd(ipntr(1)),1,0.0d0,workd(ipntr(2)),1)
|
||||
c
|
||||
c %-----------------------------------------%
|
||||
c | L O O P B A C K to call DNAUPD again. |
|
||||
c %-----------------------------------------%
|
||||
c
|
||||
go to 10
|
||||
c
|
||||
else if (ido .eq. 2) then
|
||||
|
||||
call dcopy (n, workd(ipntr(1)), 1, workd(ipntr(2)), 1)
|
||||
|
||||
go to 10
|
||||
|
||||
end if
|
||||
c
|
||||
c %----------------------------------------%
|
||||
c | Either we have convergence or there is |
|
||||
c | an error. |
|
||||
c %----------------------------------------%
|
||||
c
|
||||
if ( info .lt. 0 ) then
|
||||
c
|
||||
c %--------------------------%
|
||||
c | Error message, check the |
|
||||
c | documentation in DNAUPD. |
|
||||
c %--------------------------%
|
||||
c
|
||||
print *, ' '
|
||||
print *, ' Error with _naupd, info = ', info
|
||||
print *, ' Check the documentation of _naupd'
|
||||
print *, ' '
|
||||
stop 1
|
||||
c
|
||||
else
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | No fatal errors occurred. |
|
||||
c | Post-Process using DNEUPD. |
|
||||
c | |
|
||||
c | Computed eigenvalues may be extracted. |
|
||||
c | |
|
||||
c | Eigenvectors may also be computed now if |
|
||||
c | desired. (indicated by rvec = .true.) |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
rvec = .true.
|
||||
c
|
||||
call dneupd ( rvec, 'A', select, d, d(1,2), v, ldv,
|
||||
& sigmar, sigmai, workev, bmat, n, which, nev, tol,
|
||||
& resid, ncv, v, ldv, iparam, ipntr, workd, workl,
|
||||
& lworkl, ierr )
|
||||
c
|
||||
c %-----------------------------------------------%
|
||||
c | The real part of the eigenvalue is returned |
|
||||
c | in the first column of the two dimensional |
|
||||
c | array D, and the imaginary part is returned |
|
||||
c | in the second column of D. The corresponding |
|
||||
c | eigenvectors are returned in the first NEV |
|
||||
c | columns of the two dimensional array V if |
|
||||
c | requested. Otherwise, an orthogonal basis |
|
||||
c | for the invariant subspace corresponding to |
|
||||
c | the eigenvalues in D is returned in V. |
|
||||
c %-----------------------------------------------%
|
||||
c
|
||||
if ( ierr .ne. 0) then
|
||||
c
|
||||
c %------------------------------------%
|
||||
c | Error condition: |
|
||||
c | Check the documentation of DNEUPD. |
|
||||
c %------------------------------------%
|
||||
c
|
||||
print *, ' '
|
||||
print *, ' Error with _neupd, info = ', ierr
|
||||
print *, ' Check the documentation of _neupd. '
|
||||
print *, ' '
|
||||
c
|
||||
else
|
||||
c
|
||||
first = .true.
|
||||
nconv = iparam(5)
|
||||
do 20 j=1, nconv
|
||||
c
|
||||
c %---------------------------%
|
||||
c | Compute the residual norm |
|
||||
c | |
|
||||
c | || A*x - lambda*x || |
|
||||
c | |
|
||||
c | for the NCONV accurately |
|
||||
c | computed eigenvalues and |
|
||||
c | eigenvectors. (iparam(5) |
|
||||
c | indicates how many are |
|
||||
c | accurate to the requested |
|
||||
c | tolerance) |
|
||||
c %---------------------------%
|
||||
c
|
||||
if (d(j,2) .eq. zero) then
|
||||
c
|
||||
c %--------------------%
|
||||
c | Ritz value is real |
|
||||
c %--------------------%
|
||||
c
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& v(1,j),1,0.0d0,ax,1)
|
||||
call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
|
||||
d(j,3) = dnrm2(n, ax, 1)
|
||||
d(j,3) = d(j,3) / abs(d(j,1))
|
||||
c
|
||||
else if (first) then
|
||||
c
|
||||
c %------------------------%
|
||||
c | Ritz value is complex. |
|
||||
c | Residual of one Ritz |
|
||||
c | value of the conjugate |
|
||||
c | pair is computed. |
|
||||
c %------------------------%
|
||||
c
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& v(1,j),1,0.0d0,ax,1)
|
||||
call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
|
||||
call daxpy(n, d(j,2), v(1,j+1), 1, ax, 1)
|
||||
d(j,3) = dnrm2(n, ax, 1)
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& v(1,j+1),1,0.0d0,ax,1)
|
||||
call daxpy(n, -d(j,2), v(1,j), 1, ax, 1)
|
||||
call daxpy(n, -d(j,1), v(1,j+1), 1, ax, 1)
|
||||
d(j,3) = dlapy2( d(j,3), dnrm2(n, ax, 1) )
|
||||
d(j,3) = d(j,3) / dlapy2(d(j,1),d(j,2))
|
||||
d(j+1,3) = d(j,3)
|
||||
first = .false.
|
||||
else
|
||||
first = .true.
|
||||
end if
|
||||
c
|
||||
20 continue
|
||||
c
|
||||
c %-----------------------------%
|
||||
c | Display computed residuals. |
|
||||
c %-----------------------------%
|
||||
c
|
||||
call dmout(6, nconv, 3, d, maxncv, -6,
|
||||
& 'Ritz values (Real,Imag) and relative residuals')
|
||||
end if
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | Print additional convergence information. |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
if ( info .eq. 1) then
|
||||
print *, ' '
|
||||
print *, ' Maximum number of iterations reached.'
|
||||
print *, ' '
|
||||
else if ( info .eq. 3) then
|
||||
print *, ' '
|
||||
print *, ' No shifts could be applied during implicit',
|
||||
& ' Arnoldi update, try increasing NCV.'
|
||||
print *, ' '
|
||||
end if
|
||||
c
|
||||
print *, ' '
|
||||
print *, ' _NDRV1 '
|
||||
print *, ' ====== '
|
||||
print *, ' '
|
||||
print *, ' Size of the matrix is ', n
|
||||
print *, ' The number of Ritz values requested is ', nev
|
||||
print *, ' The number of Arnoldi vectors generated',
|
||||
& ' (NCV) is ', ncv
|
||||
print *, ' What portion of the spectrum: ', which
|
||||
print *, ' The number of converged Ritz values is ',
|
||||
& nconv
|
||||
print *, ' The number of Implicit Arnoldi update',
|
||||
& ' iterations taken is ', iparam(3)
|
||||
print *, ' The number of OP*x is ', iparam(9)
|
||||
print *, ' The convergence criterion is ', tol
|
||||
print *, ' '
|
||||
c
|
||||
end if
|
||||
c
|
||||
c %---------------------------%
|
||||
c | Done with program dndrv1. |
|
||||
c %---------------------------%
|
||||
c
|
||||
9000 continue
|
||||
c
|
||||
end
|
||||
c
|
||||
@@ -0,0 +1,414 @@
|
||||
program bug_142_gen
|
||||
c
|
||||
c Avoid taking the initial vector in the range of OP after a restart
|
||||
c (generalized case)
|
||||
c
|
||||
c Example program to illustrate the idea of reverse communication
|
||||
c for a standard nonsymmetric eigenvalue problem.
|
||||
c
|
||||
c We implement example one of ex-nonsym.doc in DOCUMENTS directory
|
||||
c
|
||||
c\Example-1
|
||||
c ... Suppose we want to solve A*x = lambda*x in regular mode,
|
||||
c where A is obtained from the standard central difference
|
||||
c discretization of the convection-diffusion operator
|
||||
c (Laplacian u) + rho*(du / dx)
|
||||
c on the unit square [0,1]x[0,1] with zero Dirichlet boundary
|
||||
c condition.
|
||||
c
|
||||
c ... OP = A and B = I.
|
||||
c
|
||||
c ... Assume "call av (nx,x,y)" computes y = A*x.c
|
||||
c
|
||||
c ... Use mode 1 of DNAUPD.
|
||||
c
|
||||
c\BeginLib
|
||||
c
|
||||
c\Routines called:
|
||||
c dnaupd ARPACK reverse communication interface routine.
|
||||
c dneupd ARPACK routine that returns Ritz values and (optionally)
|
||||
c Ritz vectors.
|
||||
c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
|
||||
c daxpy Level 1 BLAS that computes y <- alpha*x+y.
|
||||
c dnrm2 Level 1 BLAS that computes the norm of a vector.
|
||||
c av Matrix vector multiplication routine that computes A*x.
|
||||
c tv Matrix vector multiplication routine that computes T*x,
|
||||
c where T is a tridiagonal matrix. It is used in routine
|
||||
c av.
|
||||
c
|
||||
c\Author
|
||||
c Richard Lehoucq
|
||||
c Danny Sorensen
|
||||
c Chao Yang
|
||||
c Dept. of Computational &
|
||||
c Applied Mathematics
|
||||
c Rice University
|
||||
c Houston, Texas
|
||||
c
|
||||
c\SCCS Information: @(#)
|
||||
c FILE: ndrv1.F SID: 2.5 DATE OF SID: 10/17/00 RELEASE: 2
|
||||
c
|
||||
c\Remarks
|
||||
c 1. None
|
||||
c
|
||||
c\EndLib
|
||||
c---------------------------------------------------------------------------
|
||||
c
|
||||
c %-----------------------------%
|
||||
c | Define maximum dimensions |
|
||||
c | for all arrays. |
|
||||
c | MAXN: Maximum dimension |
|
||||
c | of the A allowed. |
|
||||
c | MAXNEV: Maximum NEV allowed |
|
||||
c | MAXNCV: Maximum NCV allowed |
|
||||
c %-----------------------------%
|
||||
c
|
||||
integer maxn, maxnev, maxncv, ldv
|
||||
parameter (maxn=256, maxnev=12, maxncv=30, ldv=maxn)
|
||||
c
|
||||
c %--------------%
|
||||
c | Local Arrays |
|
||||
c %--------------%
|
||||
c
|
||||
integer iparam(11), ipntr(14)
|
||||
logical select(maxncv)
|
||||
Double precision
|
||||
& ax(maxn), d(maxncv,3), resid(maxn),
|
||||
& v(ldv,maxncv), workd(3*maxn),
|
||||
& workev(3*maxncv),
|
||||
& workl(3*maxncv*maxncv+6*maxncv), a(maxn, maxn)
|
||||
c
|
||||
c %---------------%
|
||||
c | Local Scalars |
|
||||
c %---------------%
|
||||
c
|
||||
character bmat*1, which*2
|
||||
integer ido, n, nx, nev, ncv, lworkl, info, j,
|
||||
& ierr, nconv, maxitr, ishfts, mode
|
||||
Double precision
|
||||
& tol, sigmar, sigmai
|
||||
logical first, rvec
|
||||
c
|
||||
c %------------%
|
||||
c | Parameters |
|
||||
c %------------%
|
||||
c
|
||||
Double precision
|
||||
& zero
|
||||
parameter (zero = 0.0D+0)
|
||||
c
|
||||
c %-----------------------------%
|
||||
c | BLAS & LAPACK routines used |
|
||||
c %-----------------------------%
|
||||
c
|
||||
Double precision
|
||||
& dlapy2, dnrm2
|
||||
external dlapy2, dnrm2, daxpy
|
||||
c
|
||||
c %--------------------%
|
||||
c | Intrinsic function |
|
||||
c %--------------------%
|
||||
c
|
||||
intrinsic abs
|
||||
c
|
||||
c %-----------------------%
|
||||
c | Executable Statements |
|
||||
c %-----------------------%
|
||||
c
|
||||
c %--------------------------------------------------%
|
||||
c | The number NX is the number of interior points |
|
||||
c | in the discretization of the 2-dimensional |
|
||||
c | convection-diffusion operator on the unit |
|
||||
c | square with zero Dirichlet boundary condition. |
|
||||
c | The number N(=NX*NX) is the dimension of the |
|
||||
c | matrix. A standard eigenvalue problem is |
|
||||
c | solved (BMAT = 'I'). NEV is the number of |
|
||||
c | eigenvalues to be approximated. The user can |
|
||||
c | modify NX, NEV, NCV, WHICH to solve problems of |
|
||||
c | different sizes, and to get different parts of |
|
||||
c | the spectrum. However, The following |
|
||||
c | conditions must be satisfied: |
|
||||
c | N <= MAXN |
|
||||
c | NEV <= MAXNEV |
|
||||
c | NEV + 2 <= NCV <= MAXNCV |
|
||||
c %--------------------------------------------------%
|
||||
c
|
||||
nx = 10
|
||||
n = 11
|
||||
nev = 1
|
||||
ncv = 11
|
||||
do i = 1,n
|
||||
do j = 1,n
|
||||
a(i,j) = 0.15d0/11
|
||||
end do
|
||||
end do
|
||||
do j = 2,n
|
||||
a(1,j) = a(1,j) + 0.85d0
|
||||
end do
|
||||
do i = 2,n
|
||||
a(i,1) = (1-a(1,1))/10
|
||||
end do
|
||||
|
||||
if ( n .gt. maxn ) then
|
||||
print *, ' ERROR with _NDRV1: N is greater than MAXN '
|
||||
go to 9000
|
||||
else if ( nev .gt. maxnev ) then
|
||||
print *, ' ERROR with _NDRV1: NEV is greater than MAXNEV '
|
||||
go to 9000
|
||||
else if ( ncv .gt. maxncv ) then
|
||||
print *, ' ERROR with _NDRV1: NCV is greater than MAXNCV '
|
||||
go to 9000
|
||||
end if
|
||||
bmat = 'G'
|
||||
which = 'LM'
|
||||
c
|
||||
c %-----------------------------------------------------%
|
||||
c | The work array WORKL is used in DNAUPD as |
|
||||
c | workspace. Its dimension LWORKL is set as |
|
||||
c | illustrated below. The parameter TOL determines |
|
||||
c | the stopping criterion. If TOL<=0, machine |
|
||||
c | precision is used. The variable IDO is used for |
|
||||
c | reverse communication, and is initially set to 0. |
|
||||
c | Setting INFO=0 indicates that a random vector is |
|
||||
c | generated in DNAUPD to start the Arnoldi iteration. |
|
||||
c %-----------------------------------------------------%
|
||||
c
|
||||
lworkl = 3*ncv**2+6*ncv
|
||||
tol = zero
|
||||
ido = 0
|
||||
info = 0
|
||||
c
|
||||
c %---------------------------------------------------%
|
||||
c | This program uses exact shifts with respect to |
|
||||
c | the current Hessenberg matrix (IPARAM(1) = 1). |
|
||||
c | IPARAM(3) specifies the maximum number of Arnoldi |
|
||||
c | iterations allowed. Mode 1 of DNAUPD is used |
|
||||
c | (IPARAM(7) = 1). All these options can be changed |
|
||||
c | by the user. For details see the documentation in |
|
||||
c | DNAUPD. |
|
||||
c %---------------------------------------------------%
|
||||
c
|
||||
ishfts = 1
|
||||
maxitr = 300
|
||||
mode = 2
|
||||
c
|
||||
iparam(1) = ishfts
|
||||
iparam(3) = maxitr
|
||||
iparam(7) = mode
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | M A I N L O O P (Reverse communication) |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
10 continue
|
||||
c
|
||||
c %---------------------------------------------%
|
||||
c | Repeatedly call the routine DNAUPD and take |
|
||||
c | actions indicated by parameter IDO until |
|
||||
c | either convergence is indicated or maxitr |
|
||||
c | has been exceeded. |
|
||||
c %---------------------------------------------%
|
||||
c
|
||||
call dnaupd ( ido, bmat, n, which, nev, tol, resid,
|
||||
& ncv, v, ldv, iparam, ipntr, workd, workl, lworkl,
|
||||
& info )
|
||||
c
|
||||
if (ido .eq. -1 .or. ido .eq. 1) then
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | Perform matrix vector multiplication |
|
||||
c | y <--- OP*x |
|
||||
c | The user should supply his/her own |
|
||||
c | matrix vector multiplication routine here |
|
||||
c | that takes workd(ipntr(1)) as the input |
|
||||
c | vector, and return the matrix vector |
|
||||
c | product to workd(ipntr(2)). |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& workd(ipntr(1)),1,0.0d0,workd(ipntr(2)),1)
|
||||
c
|
||||
c %-----------------------------------------%
|
||||
c | L O O P B A C K to call DNAUPD again. |
|
||||
c %-----------------------------------------%
|
||||
c
|
||||
go to 10
|
||||
c
|
||||
else if (ido .eq. 2) then
|
||||
|
||||
call dcopy (n, workd(ipntr(1)), 1, workd(ipntr(2)), 1)
|
||||
|
||||
go to 10
|
||||
|
||||
end if
|
||||
c
|
||||
c %----------------------------------------%
|
||||
c | Either we have convergence or there is |
|
||||
c | an error. |
|
||||
c %----------------------------------------%
|
||||
c
|
||||
if ( info .lt. 0 ) then
|
||||
c
|
||||
c %--------------------------%
|
||||
c | Error message, check the |
|
||||
c | documentation in DNAUPD. |
|
||||
c %--------------------------%
|
||||
c
|
||||
print *, ' '
|
||||
print *, ' Error with _naupd, info = ', info
|
||||
print *, ' Check the documentation of _naupd'
|
||||
print *, ' '
|
||||
stop 1
|
||||
c
|
||||
else
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | No fatal errors occurred. |
|
||||
c | Post-Process using DNEUPD. |
|
||||
c | |
|
||||
c | Computed eigenvalues may be extracted. |
|
||||
c | |
|
||||
c | Eigenvectors may also be computed now if |
|
||||
c | desired. (indicated by rvec = .true.) |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
rvec = .true.
|
||||
c
|
||||
call dneupd ( rvec, 'A', select, d, d(1,2), v, ldv,
|
||||
& sigmar, sigmai, workev, bmat, n, which, nev, tol,
|
||||
& resid, ncv, v, ldv, iparam, ipntr, workd, workl,
|
||||
& lworkl, ierr )
|
||||
c
|
||||
c %-----------------------------------------------%
|
||||
c | The real part of the eigenvalue is returned |
|
||||
c | in the first column of the two dimensional |
|
||||
c | array D, and the imaginary part is returned |
|
||||
c | in the second column of D. The corresponding |
|
||||
c | eigenvectors are returned in the first NEV |
|
||||
c | columns of the two dimensional array V if |
|
||||
c | requested. Otherwise, an orthogonal basis |
|
||||
c | for the invariant subspace corresponding to |
|
||||
c | the eigenvalues in D is returned in V. |
|
||||
c %-----------------------------------------------%
|
||||
c
|
||||
if ( ierr .ne. 0) then
|
||||
c
|
||||
c %------------------------------------%
|
||||
c | Error condition: |
|
||||
c | Check the documentation of DNEUPD. |
|
||||
c %------------------------------------%
|
||||
c
|
||||
print *, ' '
|
||||
print *, ' Error with _neupd, info = ', ierr
|
||||
print *, ' Check the documentation of _neupd. '
|
||||
print *, ' '
|
||||
c
|
||||
else
|
||||
c
|
||||
first = .true.
|
||||
nconv = iparam(5)
|
||||
do 20 j=1, nconv
|
||||
c
|
||||
c %---------------------------%
|
||||
c | Compute the residual norm |
|
||||
c | |
|
||||
c | || A*x - lambda*x || |
|
||||
c | |
|
||||
c | for the NCONV accurately |
|
||||
c | computed eigenvalues and |
|
||||
c | eigenvectors. (iparam(5) |
|
||||
c | indicates how many are |
|
||||
c | accurate to the requested |
|
||||
c | tolerance) |
|
||||
c %---------------------------%
|
||||
c
|
||||
if (d(j,2) .eq. zero) then
|
||||
c
|
||||
c %--------------------%
|
||||
c | Ritz value is real |
|
||||
c %--------------------%
|
||||
c
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& v(1,j),1,0.0d0,ax,1)
|
||||
call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
|
||||
d(j,3) = dnrm2(n, ax, 1)
|
||||
d(j,3) = d(j,3) / abs(d(j,1))
|
||||
c
|
||||
else if (first) then
|
||||
c
|
||||
c %------------------------%
|
||||
c | Ritz value is complex. |
|
||||
c | Residual of one Ritz |
|
||||
c | value of the conjugate |
|
||||
c | pair is computed. |
|
||||
c %------------------------%
|
||||
c
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& v(1,j),1,0.0d0,ax,1)
|
||||
call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
|
||||
call daxpy(n, d(j,2), v(1,j+1), 1, ax, 1)
|
||||
d(j,3) = dnrm2(n, ax, 1)
|
||||
call dgemv('N',n,n,1.0d0,a,maxn,
|
||||
& v(1,j+1),1,0.0d0,ax,1)
|
||||
call daxpy(n, -d(j,2), v(1,j), 1, ax, 1)
|
||||
call daxpy(n, -d(j,1), v(1,j+1), 1, ax, 1)
|
||||
d(j,3) = dlapy2( d(j,3), dnrm2(n, ax, 1) )
|
||||
d(j,3) = d(j,3) / dlapy2(d(j,1),d(j,2))
|
||||
d(j+1,3) = d(j,3)
|
||||
first = .false.
|
||||
else
|
||||
first = .true.
|
||||
end if
|
||||
c
|
||||
20 continue
|
||||
c
|
||||
c %-----------------------------%
|
||||
c | Display computed residuals. |
|
||||
c %-----------------------------%
|
||||
c
|
||||
call dmout(6, nconv, 3, d, maxncv, -6,
|
||||
& 'Ritz values (Real,Imag) and relative residuals')
|
||||
end if
|
||||
c
|
||||
c %-------------------------------------------%
|
||||
c | Print additional convergence information. |
|
||||
c %-------------------------------------------%
|
||||
c
|
||||
if ( info .eq. 1) then
|
||||
print *, ' '
|
||||
print *, ' Maximum number of iterations reached.'
|
||||
print *, ' '
|
||||
else if ( info .eq. 3) then
|
||||
print *, ' '
|
||||
print *, ' No shifts could be applied during implicit',
|
||||
& ' Arnoldi update, try increasing NCV.'
|
||||
print *, ' '
|
||||
end if
|
||||
c
|
||||
print *, ' '
|
||||
print *, ' _NDRV1 '
|
||||
print *, ' ====== '
|
||||
print *, ' '
|
||||
print *, ' Size of the matrix is ', n
|
||||
print *, ' The number of Ritz values requested is ', nev
|
||||
print *, ' The number of Arnoldi vectors generated',
|
||||
& ' (NCV) is ', ncv
|
||||
print *, ' What portion of the spectrum: ', which
|
||||
print *, ' The number of converged Ritz values is ',
|
||||
& nconv
|
||||
print *, ' The number of Implicit Arnoldi update',
|
||||
& ' iterations taken is ', iparam(3)
|
||||
print *, ' The number of OP*x is ', iparam(9)
|
||||
print *, ' The convergence criterion is ', tol
|
||||
print *, ' '
|
||||
c
|
||||
end if
|
||||
c
|
||||
c %---------------------------%
|
||||
c | Done with program dndrv1. |
|
||||
c %---------------------------%
|
||||
c
|
||||
9000 continue
|
||||
c
|
||||
end
|
||||
c
|
||||
@@ -116,7 +116,7 @@ int zn() {
|
||||
int select[ncv];
|
||||
double _Complex z[(N+1)*(nev+1)];
|
||||
BLASINT ldz = N+1;
|
||||
double sigma=0;
|
||||
double _Complex sigma=0. + I*0.;
|
||||
int k;
|
||||
for (k=0; k < 3*N; ++k )
|
||||
workd[k] = 0;
|
||||
|
||||
@@ -120,7 +120,7 @@ void complex_symmetric_runner() {
|
||||
BLASINT const lworkl = 3 * (ncv * ncv) + 6 * ncv;
|
||||
|
||||
Real const tol = 0.0;
|
||||
Real const sigma = 0.0;
|
||||
std::complex<Real> const sigma(0.0, 0.0);
|
||||
|
||||
bool const rvec = true;
|
||||
|
||||
|
||||
@@ -24,6 +24,9 @@ if test x"$INTERFACE64" == x"1"; then
|
||||
AX_CHECK_COMPILE_FLAG(-fdefault-integer-8, FFLAGS="$FFLAGS -fdefault-integer-8",
|
||||
AX_CHECK_COMPILE_FLAG(-i8, FFLAGS="$FFLAGS -i8",
|
||||
AC_MSG_WARN([configure does not know how to make your Fortran compiler use 64-bit integers: set it manually via FFLAGS.])))
|
||||
AX_CHECK_COMPILE_FLAG(-fdefault-integer-8, FCFLAGS="$FCFLAGS -fdefault-integer-8",
|
||||
AX_CHECK_COMPILE_FLAG(-i8, FCFLAGS="$FCFLAGS -i8",
|
||||
AC_MSG_WARN([configure does not know how to make your Fortran compiler use 64-bit integers: set it manually via FCFLAGS.])))
|
||||
CFLAGS="$CFLAGS -DBLASINT=int64_t"
|
||||
fi
|
||||
|
||||
|
||||
@@ -0,0 +1,13 @@
|
||||
c Initialisation of the debug common block to "no debug".
|
||||
block data dbgini
|
||||
common /debug/ logfil, ndigit, mgetv0,
|
||||
& msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,
|
||||
& mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,
|
||||
& mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
|
||||
data logfil, ndigit, mgetv0,
|
||||
& msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,
|
||||
& mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,
|
||||
& mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
|
||||
& / 6, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
||||
& 0, 0, 0, 0, 0, 0 /
|
||||
end
|
||||
@@ -1,12 +0,0 @@
|
||||
! Initialisation of the debug common block to "no debug".
|
||||
block data debug_init
|
||||
common /debug/ logfil, ndigit, mgetv0, &
|
||||
msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
|
||||
mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
|
||||
mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
|
||||
data logfil, ndigit, mgetv0, &
|
||||
msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
|
||||
mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
|
||||
mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd &
|
||||
/ 6, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /
|
||||
end block data debug_init
|
||||
@@ -0,0 +1,19 @@
|
||||
c Initialisation of the stat common block.
|
||||
block data staini
|
||||
common /timing/
|
||||
& nopx, nbx, nrorth, nitref, nrstrt,
|
||||
& tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv,
|
||||
& tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv,
|
||||
& tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv,
|
||||
& tmvopx, tmvbx, tgetv0, titref, trvec
|
||||
data nopx, nbx, nrorth, nitref, nrstrt,
|
||||
& tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv,
|
||||
& tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv,
|
||||
& tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv,
|
||||
& tmvopx, tmvbx, tgetv0, titref, trvec
|
||||
& / 0, 0, 0, 0, 0,
|
||||
& 0., 0., 0., 0., 0., 0., 0.,
|
||||
& 0., 0., 0., 0., 0., 0., 0.,
|
||||
& 0., 0., 0., 0., 0., 0., 0.,
|
||||
& 0., 0., 0., 0., 0. /
|
||||
end
|
||||
Reference in New Issue
Block a user