ParMooN
 All Classes Functions Variables Friends Pages
basictools.h
1 #ifndef MPIMIS_BASICTOOLS
2 #define MPIMIS_BASICTOOLS
3 
4 #include <math.h>
5 #ifdef __MAC64__
6 #include <malloc/malloc.h>
7 #else
8 #include <malloc.h>
9 #endif
10 #include <stdio.h>
11 #include <memory.h>
12 #include <stdlib.h>
13 #include <assert.h>
14 #include "exported.h"
15 
16 typedef bool boolean_function(int, int, void*, int n); // For ns1k_scf used only
17 typedef double function_1D(double);
18 typedef double function_2D(double, double);
19 
20 typedef enum {
21  EM_DEFAULT,
22  EM_TRANSPOSED,
23  EM_LOWER_TRIANGULAR,
24  EM_UPPER_TRIANGULAR
25 } EvalMode;
26 
27 typedef enum {
28  DS_POINTWISE,
29  DS_MEANVALUE,
30  DS_LINEAR_FE,
31  DS_L2_CONSTANT
32 } Discretization_Scheme;
33 
34 #define ZERO_THRESHOLD 0.00000000000001 /* If |x| < ZERO_THRESHOLD then x is considered zero */
35 #define SQRT2 1.414213562373095049
36 #define TWOPI 6.283185307179586476
37 #define FOURPI 12.566370614359172953
38 
39 /**************************************************************************
40 *********************** Functions working with vectors ********************
41 **************************************************************************/
42 
43 /**************************************************************************
44 ****************** Norms, scalar product, orthonormalisation **************
45 **************************************************************************/
46 
47 double scalar_product(double* a, double* b, int n);
48 double euclidean_norm(double* a, int n);
49 double frobenius_norm(double* A, int rows, int cols); // TODO: Make inline (?)
50 double normalize_vector(double* a, int n);
51 
52 // returns the space distance (angle) between to directions
53 double vector_distance(double* a, double* b, int n);
54 
55 void make_orthonormal_to_system(double* a, double* A, int k, int n);
56 
57 /**************************************************************************
58 *********************** Standard arithmetic operations ********************
59 **************************************************************************/
60 
61 // Initialization
62 inline double* allocate_vector(int size) {
63  return (double*)malloc(size * sizeof(double));
64 }
65 inline void free_vector(double* v) {
66  if(v != NULL) free((void*)v);
67 }
68 inline void fill_vector(double* v, int n, double value) {
69  for(int i = 0; i < n; i++) v[i] = value;
70 }
71 inline void fill_vector_rand(double* v, int n) {
72  for(int i = 0; i < n; i++) v[i] = rand();
73 }
74 inline void clear_vector(double* v, int n) {
75  fill_vector(v, n, 0);
76 }
77 
78 // Standard arithmetics
79 inline void scale_vector(double* v, int n, double c) {
80  for(int i = 0; i < n; i++) v[i] *= c;
81 }
82 inline double sum_vector(double* v, int n);
83 inline double average_vector(double* v, int n);
84 
85 // min-max functions
86 inline double max(double a, double b){return a > b ? a : b;};
87 inline double min(double a, double b){return a < b ? a : b;};
88 inline int max(int a, int b){return a > b ? a : b;};
89 inline int min(int a, int b){return a < b ? a : b;};
90 
91 double min_vector(double* v, int n);
92 double max_vector(double* v, int n);
93 double min_abs_vector(double* v, int n);
94 double max_abs_vector(double* v, int n);
95 
96 int min_vector(int* v, int n);
97 int max_vector(int* v, int n);
98 int min_abs_vector(int* v, int n);
99 int max_abs_vector(int* v, int n);
100 
101 // Algebraic operations
102 // C += alpha * op(C) * op(B) + beta * C
103 void matrix_times_matrix(double* A, char* A_transposed, int A_rows, int A_cols, double* B, char* B_transposed, int B_rows, int B_cols, double* C, double alpha, double beta);
104 
105 // c = alpha * op(A) * b + beta * c
106 void matrix_times_vector(double* A, char* A_transposed, int A_rows, int A_cols, double* b, double* c, double alpha, double beta);
107 
108 // X = A * X
109 void triangular_matrix_times_matrix(double* A, char* uplo, double* X, int rows, int cols);
110 
111 // x = A * x
112 void triangular_matrix_times_vector(double* A, char* uplo, int n, double* x);
113 
114 // A is m x n fullmatrix, on output is itself's transposed
115 void transpose_matrix(double* A, int m, int n);
116 
117 // B is the corresponding (full) block of matrix A
118 void get_matrix_block(double* A, int rows, int cols, int block_row_start, int block_col_start, int block_rows, int block_cols, double* B);
119 
120 // A is n x n square fullmatrix, on output is itself's inverse
121 // return 0, if everything is ok, otherwise error number
122 int invert_matrix(double* A, int n);
123 
124 // Printing
125 void print_vector(double* v, int n);
126 void print_matrix(double* A, int m, int n);
127 void fprint_vector(char* filename, double* v, int n);
128 void fscan_vector(char* filename, double* v, int n);
129 void fprint_matrix(char* filename, double* A, int rows, int cols);
130 
131 // Fills the entries of the matrix via given 2D function and mesh
132 void fill_matrix(double* A, double (*)(int, int, void*), int rows, int cols, void* data);
133 void fill_matrix(double* A, function_2D* k, int rows, int cols, double* grdpoints);
134 
135 /**************************************************************************
136 **************************** Fast calculations ****************************
137 **************************************************************************/
138 
139 // return x^y
140 double ipow(double x, int y);
141 
142 // return log_a(x) (a is the base)
143 double log(double x, double a);
144 
145 // return x^y
146 int iipow(int x, int y);
147 
148 double sgn(double x);
149 double sqr(double x);
150 double cube(double x);
151 double cubert(double x);
152 int factorial(int n);
153 
154 // Fill data
155 void fill_hrefined_grid(double a, double b, int cLevels, int n, double* grid);
156 void fill_uniform_grid(double a, double b, int n, double* grid);
157 
158 void fill_exponential_data(function_1D f, int cLevels, int n, double* d);
159 void fill_uniform_data(function_1D f, int n, double* d);
160 
161 void transform_grid_m2l(double& xmin, double& xmax, int n, double* x);
162 void transform_grid_l2m(double& xmin, double& xmax, int n, double* x);
163 
164 void fill_vector_from_function(function_1D f, double* grid, int n, double* v, Discretization_Scheme ds = DS_POINTWISE);
165 
166 int get_hrefined_index(double h, int n, double x);
167 int get_uniform_index(double h, double x);
168 int get_index(double* x, int n, double x0);
169 
170 /**************************************************************************
171 ****************************** Testing tools ******************************
172 **************************************************************************/
173 void generate_orthonormal_matrix(double* Q, int k, int n);
174 void generate_vector(double* v, int n);
175 void generate_normalized_vector(double* v, int n);
176 
177 /*
178 Solves a feneral tridiafonal system of equations usinf LU
179 factorisation.
180 
181 PARAMETERS
182 diaf - array of the diafonal elements
183 ldiaf - array of the lower diafonal elements
184 udiaf - array of the upper diafonal elements
185 b - rifht-hand side of the equation
186 x - (output) solution vector
187 */
188 int solve_tridiagonal_system(double* diaf, double* ldiaf, double* udiaf, double* b, double* x, int n);
189 
190 inline int iround(double x) { return (int)(x + 0.5); }
191 
192 #endif // MPIMIS_BASICTOOLS
Generated on Tue May 10 2016 09:08:01 for ParMooN by   doxygen 1.8.5