3c97f7648072a0465c52e036944210dedb8d9810
[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 // dgemm performs one of the matrix-matrix operations
16 extern int C2F(dgemm)(char* transa, char* transb, int* m, int* n, int* k, double* alpha, double* A, int* lda, double* B, int* ldb, double* beta, double* C, int* ldc); 
17
18 // dgetrf computes an LU factorization of a general M by N matrix A (double) using partial pivoting with row interchanges
19 extern int C2F(dgetrf)(int* m, int* n, double* A, int* lda, int* ipiv, int* info); 
20
21 // dgetri computes the inverse of a matrix using the LU factorization computed by dgetrf
22 extern int C2F(dgetri)(int* n, double* A, int* lda, int* ipiv, double* work, int* lworkl, int* info); 
23
24 // zgemm performs one of the matrix-matrix operations
25 extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k, doublecomplex* alpha, doublecomplex* A, int* lda,
26         doublecomplex* B, int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc);
27
28 // zgetrd computes an LU factorization of a general M by N matrix A (complex*16) using partial pivoting with row interchanges
29 extern int C2F(zgetrf)(int* m, int* n, doublecomplex* A, int* lda, int* ipiv, int* info); 
30
31 // zgetri computes the inverse of a matrix using the LU factorization computed by zgetrf
32 extern int C2F(zgetri)(int* n, doublecomplex* A, int* lda, int* ipiv, doublecomplex* work, int* lworkl, int* info);
33
34 static double alpha = 1.;
35 static double beta = 0.;
36
37 static doublecomplex alphac = {1.,0.};
38 static doublecomplex betac = {0.,0.};
39
40 void RtimesRprime(double* result, double* R, double* Rprime, int N)
41 {
42         C2F(dgemm)("n", "n", &N, &N, &N, &alpha, R, &N, Rprime, &N, &beta, result, &N);
43 }
44
45 void invR_times_A_times_invRprime(double* result, double* R, double* A, double* Rprime, int N)
46 {
47         double* invRxA = NULL;
48         int* IPVT = NULL;
49         double* work = NULL;
50         double* invR = NULL;
51         double* invRprime = NULL;
52         int INFO_LU = 0;
53
54         invRxA = (double*)malloc(N * N * sizeof(double));
55         invR = (double*)malloc(N * N * sizeof(double));
56         invRprime = (double*)malloc(N * N * sizeof(double));
57         work = (double*)malloc(N * N * sizeof(double));
58
59
60         IPVT = (int*) malloc(N * sizeof(int));
61         memset(IPVT, 0, N * sizeof(int));
62
63         // inv(R)
64         memcpy(invR, R, N * N * sizeof(double));   // copy R to invR
65         C2F(dgetrf)(&N, &N, invR ,&N, IPVT, &INFO_LU);  // LU decomposition
66
67         memset(work, 0, N * N * sizeof(double));
68         C2F(dgetri)(&N, invR, &N, IPVT, work, &N, &INFO_LU);  // Obtain inverse matrix R
69
70         // inv(Rprime)
71         memset(IPVT, 0, N * sizeof(int));
72         memcpy(invRprime, Rprime, N * N * sizeof(double));      
73         C2F(dgetrf)(&N, &N, invRprime, &N, IPVT, &INFO_LU);     // LU decomposition
74
75         memset(work, 0, N * N * sizeof(double));
76         C2F(dgetri)(&N, invRprime, &N, IPVT, work, &N,&INFO_LU);        // Obtain inverse matrix Rprime
77
78         C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invR, &N, A, &N, &beta, invRxA, &N);
79
80         C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invRxA, &N, invRprime, &N, &beta, result, &N);
81
82         free(invRxA);
83         free(IPVT);
84         free(work);
85         free(invR);
86         free(invRprime);
87 }
88
89 void invU_times_invL_times_E(double* result, double* U, double* L, double* E, int N)
90 {
91         double* invUxinvL = NULL;
92         int* IPVT = NULL;
93         double* work = NULL;
94         double* invU = NULL;
95         double* invL = NULL;
96         int INFO_LU = 0;
97
98         invUxinvL = (double*)malloc(N * N * sizeof(double));
99         invU = (double*)malloc(N * N * sizeof(double));
100         invL = (double*)malloc(N * N* sizeof(double));
101         work = (double*)malloc(N * N * sizeof(double));
102
103         IPVT = (int*) malloc(N * sizeof(int));
104         memset(IPVT, 0, N * sizeof(int));
105         
106         // inv L -> L obtained with LU decomposition
107         memcpy(invL, L, N * N * sizeof(double));
108         C2F(dgetrf)(&N, &N, invL, &N, IPVT, &INFO_LU); // LU decomposition
109
110         memset(work, 0, N * N * sizeof(double));
111         C2F(dgetri)(&N, invL, &N, IPVT, work, &N, &INFO_LU);  // inv(L)
112
113         // inv U -> U obtained with LU decomposition
114         memcpy(invU, U, N * N * sizeof(double));
115         memset(IPVT, 0, N*sizeof(int));
116         C2F(dgetrf)(&N, &N, invU, &N, IPVT, &INFO_LU); // LU decomposition
117
118         memset(work, 0, N * N * sizeof(double));
119         C2F(dgetri)(&N, invU, &N, IPVT, work, &N, &INFO_LU); // inv(U)
120
121         C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invU, &N, invL, &N, &beta, invUxinvL, &N);
122
123         C2F(dgemm)("n", "n", &N, &N, &N, &alpha, invUxinvL, &N, E, &N, &beta, result, &N);
124
125         free(invUxinvL);
126         free(IPVT);
127         free(work);
128         free(invU);
129         free(invL);
130 }
131
132
133 void RCtimesRCprime(doublecomplex* result, doublecomplex* R, doublecomplex* Rprime, int N)
134 {
135         C2F(zgemm)("n", "n", &N, &N, &N, &alphac, R, &N, Rprime, &N, &betac, result, &N);
136 }
137
138 void invRC_times_AC_times_invRCprime(doublecomplex* result, doublecomplex* R, doublecomplex* A, doublecomplex* Rprime, int N)
139 {
140         doublecomplex* invRxA = NULL;
141         int* IPVT = NULL;
142         doublecomplex* work = NULL;
143         doublecomplex* invR = NULL;
144         doublecomplex* invRprime = NULL;
145         int INFO_LU = 0;
146
147         invRxA = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
148         invR = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
149         invRprime = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
150         work = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
151
152
153         IPVT = (int*) malloc(N * sizeof(int));
154         memset(IPVT, 0, N * sizeof(int));
155
156         // inv(R)
157         memcpy(invR, R, N * N * sizeof(doublecomplex));   // copy R to invR
158         C2F(zgetrf)(&N, &N, invR ,&N, IPVT, &INFO_LU);  // LU decomposition
159
160         memset(work, 0, N * N * sizeof(doublecomplex));
161         C2F(zgetri)(&N, invR, &N, IPVT, work, &N, &INFO_LU);  // Obtain inverse matrix R
162
163         // inv(Rprime)
164         memset(IPVT, 0, N * sizeof(int));
165         memcpy(invRprime, Rprime, N * N * sizeof(doublecomplex));       
166         C2F(zgetrf)(&N, &N, invRprime, &N, IPVT, &INFO_LU);     // LU decomposition
167
168         memset(work, 0, N * N * sizeof(doublecomplex));
169         C2F(zgetri)(&N, invRprime, &N, IPVT, work, &N,&INFO_LU);        // Obtain inverse matrix Rprime
170
171         C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invR, &N, A, &N, &betac, invRxA, &N);
172
173         C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invRxA, &N, invRprime, &N, &betac, result, &N);
174
175         free(invRxA);
176         free(IPVT);
177         free(work);
178         free(invR);
179         free(invRprime);
180 }
181
182
183 void invUC_times_invLC_times_EC(doublecomplex* result, doublecomplex* U, doublecomplex* L, doublecomplex* E, int N)
184 {
185         doublecomplex* invUxinvL = NULL;
186         int* IPVT = NULL;
187         doublecomplex* work = NULL;
188         doublecomplex* invU = NULL;
189         doublecomplex* invL = NULL;
190         int INFO_LU = 0;
191
192         invUxinvL = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
193         invU = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
194         invL = (doublecomplex*)malloc(N * N* sizeof(doublecomplex));
195         work = (doublecomplex*)malloc(N * N * sizeof(doublecomplex));
196
197         IPVT = (int*) malloc(N * sizeof(int));
198         memset(IPVT, 0, N * sizeof(int));
199         
200         // inv L -> L obtained with LU decomposition
201         memcpy(invL, L, N * N * sizeof(doublecomplex));
202         C2F(zgetrf)(&N, &N, invL, &N, IPVT, &INFO_LU); // LU decomposition
203
204         memset(work, 0, N * N * sizeof(doublecomplex));
205         C2F(zgetri)(&N, invL, &N, IPVT, work, &N, &INFO_LU);  // inv(L)
206
207         // inv U -> U obtained with LU decomposition
208         memcpy(invU, U, N * N * sizeof(doublecomplex));
209         memset(IPVT, 0, N*sizeof(int));
210         C2F(zgetrf)(&N, &N, invU, &N, IPVT, &INFO_LU); // LU decomposition
211
212         memset(work, 0, N * N * sizeof(doublecomplex));
213         C2F(zgetri)(&N, invU, &N, IPVT, work, &N, &INFO_LU); // inv(U)
214
215         C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invU, &N, invL, &N, &betac, invUxinvL, &N);
216
217         C2F(zgemm)("n", "n", &N, &N, &N, &alphac, invUxinvL, &N, E, &N, &betac, result, &N);
218
219         free(invUxinvL);
220         free(IPVT);
221         free(work);
222         free(invU);
223         free(invL);
224 }