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