Bug #12138 fixed - eigs(A,B) returned incorrect eigenvectors for dense matrices.
[scilab.git] / 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 }