Add unitary tests for arnoldi
[scilab.git] / scilab / modules / arnoldi / tests / unit_tests / dseupd.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 d       = zeros(nev, 1);
21 z       = zeros(nx, nev);
22 resid   = zeros(nx, 1); 
23 v       = zeros(nx, ncv);
24 workd   = zeros(3 * nx, 1); 
25 workl   = zeros(ncv * ncv + 8 * ncv, 1);
26
27 // Build the symmetric test matrix
28
29 A            = diag(10 * ones(nx,1));
30 A(1:$-1,2:$) = A(1:$-1,2:$) + diag(6 * ones(nx-1,1));
31 A(2:$,1:$-1) = A(2:$,1:$-1) + diag(6 * ones(nx-1,1));
32
33 tol    = 0;
34 ido    = 0;
35
36 ishfts = 1;
37 maxitr = 300;
38 mode1  = 1;
39
40 iparam(1) = ishfts;
41 iparam(3) = maxitr;
42 iparam(7) = mode1;
43
44 sigma = 0; // the real part of the shift
45 info_dsaupd = 0;
46
47 // M A I N   L O O P (Reverse communication)
48
49 while(ido <> 99)
50   // Repeatedly call the routine DSAUPD and take actions indicated by parameter IDO until
51   // either convergence is indicated or maxitr has been exceeded.
52
53   [ido, resid, v, iparam, ipntr, workd, workl, info_dsaupd] = dsaupd(ido, bmat, nx, which, nev, tol, resid, ncv, v, iparam, ipntr, workd, workl, info_dsaupd);
54   
55   if(info_dsaupd < 0)
56     printf('\nError with dsaupd, info = %d\n',info_dsaupd);
57     printf('Check the documentation of dsaupd\n\n');
58   end
59   
60   if(ido == -1 | ido == 1)
61     // Perform matrix vector multiplication 
62     workd(ipntr(2):ipntr(2) + nx - 1) = A * workd(ipntr(1):ipntr(1) + nx - 1);
63   end
64 end
65
66 // Post-Process using DSEUPD.
67 rvec    = 1;
68 howmany = 'A';
69 info_dseupd = 0;
70
71 [d, z, resid, v, iparam, ipntr, workd, workl, info_dseupd] = dseupd(rvec, howmany, _select, d, z, sigma, bmat, nx, which, nev, tol, resid, ncv, v, ...
72                                                                     iparam, ipntr, workd, workl, info_dseupd);
73
74 assert_checkalmostequal(A * z, z * diag(d), sqrt(%eps), 1.e-10);