MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
arptest.cc File Reference
#include "mbconfig.h"
#include <iostream>
#include <iomanip>
#include <vector>
#include "ac/f2c.h"
#include "ac/arpack.h"
#include "solman.h"
#include "naivewrap.h"
Include dependency graph for arptest.cc:

Go to the source code of this file.

Functions

int main (void)
 

Function Documentation

int main ( void  )

Definition at line 742 of file arptest.cc.

References NaiveMatrixHandler::iGetNumRows(), and MatrixHandler::MatVecIncMul().

743 {
744  // matrix
745  NaiveMatrixHandler mh(5);
746  mh(1, 1) = 2.;
747  mh(1, 2) = 1.;
748  mh(2, 1) = 1.;
749  mh(2, 2) = 2.;
750  mh(2, 3) = 1.;
751  mh(3, 2) = 1.;
752  mh(3, 3) = 2.;
753  mh(3, 4) = 1.;
754  mh(4, 3) = 1.;
755  mh(4, 4) = 2.;
756  mh(4, 5) = 1.;
757  mh(5, 4) = 1.;
758  mh(5, 5) = 2.;
759 
760  // shift
761  doublereal SIGMA = 1.;
762  doublereal SIGMAR = 0.;
763  doublereal SIGMAI = 0.;
764 
765  // arpack-related vars
766  integer IDO; // 0 at first iteration; then set by dnaupd
767  const char *BMAT; // 'I' for standard problem
768  integer N; // size of problem
769  const char *WHICH; // "SM" to request smallest eigenvalues
770  integer NEV; // number of eigenvalues
771  doublereal TOL; // -1 to use machine precision
772  std::vector<doublereal> RESID; // residual vector (ignored if IDO==0)
773  integer NCV; // number of vectors in subspace
774  std::vector<doublereal> V; // Schur basis
775  integer LDV; // leading dimension of V (==N!)
776  integer IPARAM[11] = { 0 };
777  integer IPNTR[14] = { 0 };
778  std::vector<doublereal> WORKD;
779  std::vector<doublereal> WORKL;
780  integer LWORKL;
781  integer INFO;
782 
783  IDO = 0;
784  BMAT = "I";
785  N = mh.iGetNumRows();
786  WHICH = "SM";
787  NEV = 2;
788  TOL = 0.;
789  RESID.resize(N, 0.);
790  NCV = 4;
791  V.resize(N*NCV, 0.);
792  LDV = N;
793  IPARAM[0] = 1;
794  IPARAM[2] = 300;
795  IPARAM[3] = 1;
796  IPARAM[6] = 1; // mode...
797  WORKD.resize(3*N, 0.);
798  LWORKL = 3*NCV*NCV + 6*NCV;
799  WORKL.resize(LWORKL, 0.);
800  INFO = 0;
801 
802  int cnt = 0;
803  do {
804  std::cout << "before:" << std::endl
805  << "IDO=" << IDO << std::endl
806  << "BMAT=" << BMAT << std::endl
807  << "N=" << N << std::endl
808  << "WHICH=" << WHICH << std::endl
809  << "NEV=" << NEV << std::endl
810  << "TOL=" << TOL << std::endl
811  << "&RESID[0]=" << &RESID[0] << std::endl
812  << "NCV=" << NCV << std::endl
813  << "&V[0]=" << &V[0] << std::endl
814  << "LDV=" << LDV << std::endl
815  << "IPARAM=" << IPARAM[0]
816  << "," << IPARAM[1]
817  << "," << IPARAM[2]
818  << "," << IPARAM[3]
819  << "," << IPARAM[4]
820  << "," << IPARAM[5]
821  << "," << IPARAM[6]
822  << "," << IPARAM[7]
823  << "," << IPARAM[8]
824  << "," << IPARAM[9]
825  << "," << IPARAM[10]
826  << std::endl
827  << "IPNTR=" << IPNTR[0]
828  << "," << IPNTR[1]
829  << "," << IPNTR[2]
830  << "," << IPNTR[3]
831  << "," << IPNTR[4]
832  << "," << IPNTR[5]
833  << "," << IPNTR[6]
834  << "," << IPNTR[7]
835  << "," << IPNTR[8]
836  << "," << IPNTR[9]
837  << "," << IPNTR[10]
838  << "," << IPNTR[11]
839  << "," << IPNTR[12]
840  << "," << IPNTR[13]
841  << std::endl
842  << "&WORKD[0]=" << &WORKD[0] << std::endl
843  << "&WORKL[0]=" << &WORKL[0] << std::endl
844  << "LWORKL=" << LWORKL << std::endl
845  << "INFO=" << INFO << std::endl
846  << std::endl;
847 
848  __FC_DECL__(dnaupd)(&IDO, &BMAT[0], &N, &WHICH[0], &NEV,
849  &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
850  &WORKD[0], &WORKL[0], &LWORKL, &INFO);
851 
852  std::cout << "after:" << std::endl
853  << "IDO=" << IDO << std::endl
854  << "BMAT=" << BMAT << std::endl
855  << "N=" << N << std::endl
856  << "WHICH=" << WHICH << std::endl
857  << "NEV=" << NEV << std::endl
858  << "TOL=" << TOL << std::endl
859  << "&RESID[0]=" << &RESID[0] << std::endl
860  << "NCV=" << NCV << std::endl
861  << "&V[0]=" << &V[0] << std::endl
862  << "LDV=" << LDV << std::endl
863  << "IPARAM=" << IPARAM[0]
864  << "," << IPARAM[1]
865  << "," << IPARAM[2]
866  << "," << IPARAM[3]
867  << "," << IPARAM[4]
868  << "," << IPARAM[5]
869  << "," << IPARAM[6]
870  << "," << IPARAM[7]
871  << "," << IPARAM[8]
872  << "," << IPARAM[9]
873  << "," << IPARAM[10]
874  << std::endl
875  << "IPNTR=" << IPNTR[0]
876  << "," << IPNTR[1]
877  << "," << IPNTR[2]
878  << "," << IPNTR[3]
879  << "," << IPNTR[4]
880  << "," << IPNTR[5]
881  << "," << IPNTR[6]
882  << "," << IPNTR[7]
883  << "," << IPNTR[8]
884  << "," << IPNTR[9]
885  << "," << IPNTR[10]
886  << "," << IPNTR[11]
887  << "," << IPNTR[12]
888  << "," << IPNTR[13]
889  << std::endl
890  << "&WORKD[0]=" << &WORKD[0] << std::endl
891  << "&WORKL[0]=" << &WORKL[0] << std::endl
892  << "LWORKL=" << LWORKL << std::endl
893  << "INFO=" << INFO << std::endl
894  << std::endl;
895 
896  std::cout << "cnt=" << cnt << ": IDO=" << IDO << ", INFO=" << INFO << std::endl;
897 
898  // compute Y = OP*X
899  MyVectorHandler X(N, &WORKD[IPNTR[0] - 1]);
900  MyVectorHandler Y(N, &WORKD[IPNTR[1] - 1]);
901 
902  Y = X;
903  Y *= -SIGMA;
904  mh.MatVecIncMul(Y, X);
905 
906  std::cout << "X:" << std::endl << X << std::endl;
907  std::cout << "Y:" << std::endl << Y << std::endl;
908 
909  cnt++;
910  } while (IDO == 1 || IDO == -1);
911 
912  if (INFO < 0) {
913  std::cerr << "error" << std::endl;
914  return 1;
915  }
916 
917  logical RVEC = true;
918  const char *HOWMNY = "A";
919  std::vector<logical> SELECT(NCV);
920  std::vector<doublereal> DR(NEV + 1);
921  std::vector<doublereal> DI(NEV + 1);
922  std::vector<doublereal> Z(N*(NEV + 1));
923  integer LDZ = N;
924  std::vector<doublereal> WORKEV(3*NCV);
925 
926  __FC_DECL__(dneupd)(&RVEC, &HOWMNY[0], &SELECT[0], &DR[0], &DI[0],
927  &Z[0], &LDZ, &SIGMAR, &SIGMAI, &WORKEV[0],
928  &BMAT[0], &N, &WHICH[0], &NEV,
929  &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
930  &WORKD[0], &WORKL[0], &LWORKL, &INFO);
931 
932  bool bCmplx = false;
933  bool bFirst;
934  for (integer nv = 0; nv < NEV; nv++) {
935  std::cout << "value " << std::setw(8) << nv << ":"
936  << std::setw(16) << DR[nv] + SIGMA
937  << (DI[nv] >= 0. ? " + " : " - ")
938  << std::setw(16) << std::abs(DI[nv]) << " i" << std::endl;
939 
940  std::cout << "vector " << std::setw(8) << nv << ":"
941  << std::endl;
942  if (!bCmplx) {
943  if (DI[nv] != 0.) {
944  bCmplx = true;
945  bFirst = true;
946  }
947  }
948 
949  if (!bCmplx) {
950  for (integer idx = 0; idx < N; idx++) {
951  std::cout
952  << std::setw(8) << idx << ":"
953  << std::setw(16) << Z[N*nv + idx]
954  << std::endl;
955  }
956 
957  } else if (bFirst) {
958  bFirst = false;
959 
960  for (integer idx = 0; idx < N; idx++) {
961  doublereal d = Z[N*(nv + 1) + idx];
962  std::cout
963  << std::setw(8) << idx << ":"
964  << std::setw(16) << Z[N*nv + idx]
965  << (d > 0. ? " + " : " - " )
966  << std::setw(16) << std::abs(d) << " i"
967  << std::endl;
968  }
969 
970  } else {
971  bCmplx = false;
972 
973  for (integer idx = 0; idx < N; idx++) {
974  doublereal d = Z[N*nv + idx];
975  std::cout
976  << std::setw(8) << idx << ":"
977  << std::setw(16) << Z[N*(nv - 1) + idx]
978  << (d > 0. ? " + " : " - " )
979  << std::setw(16) << std::abs(d) << " i"
980  << std::endl;
981  }
982  }
983  }
984 
985  return 0;
986 }
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51

Here is the call graph for this function: