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