1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
|
/**********************************************************************
ReLU_inverse.c:
ReLU_inverse.c is a subroutine to calculate the inverse of matrix
constructed by real number using LU factorization.
Log of ReLU_inverse.c:
22/Nov/2001 Released by T.Ozaki
***********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "openmx_common.h"
void ReLU_inverse(int n, double **a, double **ia)
{
static int i,j,k;
static double w;
static double *x,*y;
static double **da;
/****************************************************
allocation of arrays:
static double x[List_YOUSO[7]];
static double y[List_YOUSO[7]];
static double da[List_YOUSO[7]][List_YOUSO[7]];
****************************************************/
x = (double*)malloc(sizeof(double)*List_YOUSO[7]);
y = (double*)malloc(sizeof(double)*List_YOUSO[7]);
da = (double**)malloc(sizeof(double*)*List_YOUSO[7]);
for (i=0; i<List_YOUSO[7]; i++){
da[i] = (double*)malloc(sizeof(double)*List_YOUSO[7]);
}
/***************************************************
start calc.
****************************************************/
if (n==-1){
for (i=0; i<List_YOUSO[7]; i++){
for (j=0; j<List_YOUSO[7]; j++){
a[i][j] = 0.0;
}
}
}
else{
for (i=0; i<=n; i++){
for (j=0; j<=n; j++){
da[i][j] = a[i][j];
}
}
/****************************************************
LU factorization
****************************************************/
for (k=0; k<=n-1; k++){
w = 1.0/a[k][k];
for (i=k+1; i<=n; i++){
a[i][k] = w*a[i][k];
for (j=k+1; j<=n; j++){
a[i][j] = a[i][j] - a[i][k]*a[k][j];
}
}
}
for (k=0; k<=n; k++){
/****************************************************
Ly = b
****************************************************/
for (i=0; i<=n; i++){
if (i==k)
y[i] = 1.0;
else
y[i] = 0.0;
for (j=0; j<=i-1; j++){
y[i] = y[i] - a[i][j]*y[j];
}
}
/****************************************************
Ly = b
****************************************************/
for (i=n; 0<=i; i--){
x[i] = y[i];
for (j=n; (i+1)<=j; j--){
x[i] = x[i] - a[i][j]*x[j];
}
x[i] = x[i]/a[i][i];
ia[i][k] = x[i];
}
}
for (i=0; i<=n; i++){
for (j=0; j<=n; j++){
a[i][j] = da[i][j];
}
}
}
/****************************************************
freeing of arrays:
static double x[List_YOUSO[7]];
static double y[List_YOUSO[7]];
static double da[List_YOUSO[7]][List_YOUSO[7]];
****************************************************/
free(x);
free(y);
for (i=0; i<List_YOUSO[7]; i++){
free(da[i]);
}
free(da);
}
|