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 }