scilab/modules/arnoldi/src/c/eigs_dependencies.c

   1 /*
   2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3  * Copyright (C) 2011 - Scilab Enterprises - Adeline CARNIS
   4  *
   5  * This file must be used under the terms of the CeCILL.
   6  * This source file is licensed as described in the file COPYING, which
   7  * you should have received as part of this distribution.  The terms
   8  * are also available at
   9  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  10  *
  11  */
  12
  13 #include "eigs_dependencies.h"
  14
  15 extern int C2F(dgemv) (char* trans, int* m, int* n, double* alpha, double* A,
  16                        int* lda, double* x, int* incx, double* beta, double* y, int* incy);
  17 extern double C2F(ddot) (int *n, double* x, int* incx, double* y, int* incy);
  18
  19 void process_dneupd_data(double* DR, double* DI, double* Z, int N, int nev, double* AR,
  20                          doublecomplex* eigenvalue, doublecomplex* eigenvector,
  21                          int sigma_imaginary)
  22 {
  23     /* if sigma_imaginary there is an extra step to compute the eigenvalues
  24        as explained in the dneupd user guide */
  25
  26     double* temp1 = NULL;
  27     double* temp2 = NULL;
  28
  29     int i = 0;
  30     int j = 0;
  31
  32     double alpha = 1;
  33     double beta = 0;
  34     int iOne = 1;
  35     double real_part;
  36     double imag_part;
  37
  38     if ( sigma_imaginary )
  39     {
  40         temp1 = (double*) malloc(N * sizeof(double));
  41         temp2 = (double*) malloc(N * sizeof(double));
  42
  43         while (i < nev)
  44         {
  45             if (DI[i] == 0)
  46             {
  47                 C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * i, &iOne, &beta, temp1, &iOne);
  48                 eigenvalue[i] = (doublecomplex)
  49                 {
  50                     C2F(ddot) (&N, Z + N * i, &iOne, temp1, &iOne), 0
  51                 };
  52                 i = i + 1;
  53             }
  54             else
  55             {
  56                 C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * i, &iOne, &beta, temp1, &iOne);
  57                 C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * (i + 1), &iOne, &beta, temp2, &iOne);
  58                 real_part = C2F(ddot) (&N, Z + N * i, &iOne, temp1, &iOne) + \
  59                             C2F(ddot) (&N, Z + N * (i + 1), &iOne, temp2, &iOne);
  60                 imag_part = C2F(ddot) (&N, Z + N * i, &iOne, temp2, &iOne) - \
  61                             C2F(ddot) (&N, Z + N * (i + 1), &iOne, temp1, &iOne);
  62                 eigenvalue[i] = (doublecomplex)
  63                 {
  64                     real_part, imag_part
  65                 };
  66                 eigenvalue[i + 1] = (doublecomplex)
  67                 {
  68                     real_part, -imag_part
  69                 };
  70                 i = i + 2;
  71             }
  72         }
  73         free(temp1);
  74         free(temp2);
  75     }
  76     else
  77     {
  78         for (i = 0; i < nev + 1; i++)
  79         {
  80             eigenvalue[i] = (doublecomplex)
  81             {
  82                 DR[i], DI[i]
  83             };
  84         }
  85     }
  86
  87     if (eigenvector)
  88     {
  89         i = 0;
  90
  91         while (i < nev)
  92         {
  93             if (DI[i] != 0)
  94             {
  95                 for (j = 0; j < N; j++)
  96                 {
  97                     eigenvector[i * N + j] = (doublecomplex)
  98                     {
  99                         Z[i * N + j], Z[(i + 1) * N + j]
 100                     };
 101                     eigenvector[(i + 1) * N + j] = (doublecomplex)
 102                     {
 103                         Z[i * N + j], -Z[(i + 1) * N + j]
 104                     };
 105                 }
 106                 i = i + 2;
 107             }
 108             else
 109             {
 110                 for (j = 0; j < N; j++)
 111                 {
 112                     eigenvector[i * N + j] = (doublecomplex)
 113                     {
 114                         Z[i * N + j], 0
 115                     };
 116                 }
 117                 i = i + 1;
 118             }
 119         }
 120     }
 121 }