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