1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2012 - SE - Sylvestre Ledru
4//
5//  This file is distributed under the same license as the Scilab package.
6// =============================================================================
7
8// <-- CLI SHELL MODE -->
9
10nx    = 10;
11
12nev   = 3;
13ncv   = 6;
14bmat  = 'I';
15which = 'LM';
16
17// Local Arrays
18
19iparam  = zeros(11, 1);
20ipntr   = zeros(14, 1);
21_select = zeros(ncv, 1);
22dr       = zeros(nev + 1, 1);
23di      = zeros(nev + 1, 1);
24z       = zeros(nx, nev + 1);
25resid   = zeros(nx, 1);
26v       = zeros(nx, ncv);
27workd   = zeros(3 * nx, 1);
28workev  = zeros(3 * ncv, 1);
29workl   = zeros(3 * ncv * ncv + 6 * ncv, 1);
30
31// Build the test matrix
32
33A            = diag(10 * ones(nx, 1));
34A(1:$-1,2:$) = A(1:$-1,2:$) + diag(6 * ones(nx-1,1));
35A(2:$,1:$-1) = A(2:$,1:$-1) + diag(-6 * ones(nx-1,1));
36
37tol    = 0;
38ido    = 0;
39
40ishfts = 1;
41maxitr = 300;
42mode1  = 1;
43
44iparam(1) = ishfts;
45iparam(3) = maxitr;
46iparam(7) = mode1;
47
48sigmar = 0; // the real part of the shift
49sigmai = 0; // the imaginary part of the shift
50info_dnaupd = 0;
51
52// M A I N   L O O P (Reverse communication)
53
54while(ido <> 99)
55    // Repeatedly call the routine DNAUPD and take actions indicated by parameter IDO until
56    // either convergence is indicated or maxitr has been exceeded.
57
58    [ido, resid, v, iparam, ipntr, workd, workl, info_dnaupd] = dnaupd(ido, bmat, nx, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_dnaupd);
59
60    if(info_dnaupd < 0)
61        printf('\nError with dnaupd, info = %d\n',info_dnaupd);
62        printf('Check the documentation of dnaupd\n\n');
63    end
64
65    if(ido == -1 | ido == 1)
66        // Perform matrix vector multiplication
67        workd(ipntr(2):ipntr(2) + nx -1) = A * workd(ipntr(1):ipntr(1) + nx - 1);
68    end
69end
70
71// Post-Process using DNEUPD.
72rvec    = 1;
73howmany = 'A';
74info_dneupd = 0;
75
76[dr, di, z, resid, v, iparam, ipntr, workd, workl, info_dneupd] = dneupd(rvec, howmany, _select, dr, di, z, sigmar, sigmai, workev, ...
77bmat, nx, which, nev, tol, resid, ncv, v, ...
78iparam, ipntr, workd, workl, info_dneupd);
79
80d = complex(dr, di);
81d(nev+1) = [];
82d = diag(d);
83
84c1 = 1:2:nev + 1;
85c2 = 2:2:nev + 1;
86if(modulo(nev + 1, 2) == 1)
87    c1($) = [];
88end
89z(:,[c1, c2]) = [z(:,c1) + z(:,c2) * %i z(:,c1) - z(:,c2) * %i];
90z(:,$) = [];
91
92assert_checkalmostequal(A * z, z * d, sqrt(%eps), 1.e-10);
93