libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
matrix.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of libdspl-2.0.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include "dspl.h"
24 #include "dspl_internal.h"
25 #include "blas.h"
26 
27 
28 
29 
30 #ifdef DOXYGEN_ENGLISH
31 
32 #endif
33 #ifdef DOXYGEN_RUSSIAN
34 
89 #endif
90 int DSPL_API matrix_eig_cmplx(complex_t* a, int n, complex_t* v, int* info)
91 {
92  int err;
93  int sdim = 0;
94  int ldvs = 1;
95  int lwork = 2*n;
96  if(!a || !v)
97  return ERROR_PTR;
98 
99  if(n<1)
100  return ERROR_MATRIX_SIZE;
101 
102  complex_t *work=(complex_t*)malloc(lwork*sizeof(complex_t));
103  double *rwork = (double*)malloc(n*sizeof(double));
104 
105  zgees_("N", "N", NULL, &n, a, &n, &sdim, v, NULL, &ldvs, work, &lwork,
106  rwork, NULL, &err);
107 
108  if(err!=0)
109  {
110  if(info)
111  *info = err;
112  err = ERROR_LAPACK;
113  }
114  else
115  err = RES_OK;
116 
117  free(work);
118  free(rwork);
119  return err;
120 }
121 
122 
123 
124 #ifdef DOXYGEN_ENGLISH
125 
126 #endif
127 #ifdef DOXYGEN_RUSSIAN
128 
152 #endif
153 int DSPL_API matrix_eye(double* a, int n, int m)
154 {
155  int p, k;
156  if(!a)
157  return ERROR_PTR;
158  if (n < 1 || m < 1)
159  return ERROR_MATRIX_SIZE;
160 
161  k = 0;
162  memset(a, 0, n*m*sizeof(double));
163  for(p = 0; p < m; p++)
164  {
165  a[k] = 1.0;
166  k += n+1;
167  }
168 
169  return RES_OK;
170 }
171 
172 
173 #ifdef DOXYGEN_ENGLISH
174 
175 #endif
176 #ifdef DOXYGEN_RUSSIAN
177 
178 #endif
179 int DSPL_API matrix_eye_cmplx(complex_t* a, int n, int m)
180 {
181  int p, k;
182  if(!a)
183  return ERROR_PTR;
184  if (n < 1 || m < 1)
185  return ERROR_MATRIX_SIZE;
186 
187  k = 0;
188  memset(a, 0, n*m*sizeof(complex_t));
189  for(p = 0; p < m; p++)
190  {
191  RE(a[k]) = 1.0;
192  k += n+1;
193  }
194 
195  return RES_OK;
196 }
197 
198 
199 
200 
201 
202 #ifdef DOXYGEN_ENGLISH
203 
204 #endif
205 #ifdef DOXYGEN_RUSSIAN
206 
207 #endif
208 int DSPL_API matrix_mul(double* a, int na, int ma,
209  double* b, int nb, int mb,
210  double* c)
211 {
212 
213  double alpha = 1;
214  double beta = 0.0;
215 
216  if(!a || !b || !c)
217  return ERROR_PTR;
218  if(na < 1 || ma < 1 || nb < 1 || mb < 1 || ma != nb)
219  return ERROR_MATRIX_SIZE;
220 
221  /* BLAS DGEMM */
222  dgemm_("N", "N", &na, &mb, &ma, &alpha, a, &na, b, &nb, &beta, c, &na);
223 
224  return RES_OK;
225 }
226 
227 
228 
229 #ifdef DOXYGEN_ENGLISH
230 
231 #endif
232 #ifdef DOXYGEN_RUSSIAN
233 
234 #endif
235 int DSPL_API matrix_print(double* a, int n, int m,
236  const char* name, const char* format)
237 {
238  int p,q;
239 
240  if(!a)
241  return ERROR_PTR;
242  if(n < 1 || m < 1)
243  return ERROR_SIZE;
244 
245  printf("\n%s = [ %% size [%d x %d] type: real", name, n, m);
246 
247  for(p = 0; p < n; p++)
248  {
249  printf("\n");
250  for(q = 0; q < m; q++)
251  {
252  printf(format, a[q*n + p]);
253  if(q == m-1)
254  printf(";");
255  else
256  printf(", ");
257  }
258  }
259  printf("];\n");
260 
261  return RES_OK;
262 }
263 
264 
265 #ifdef DOXYGEN_ENGLISH
266 
267 #endif
268 #ifdef DOXYGEN_RUSSIAN
269 
270 #endif
271 int DSPL_API matrix_print_cmplx(complex_t* a, int n, int m,
272  const char* name, const char* format)
273 {
274  int p,q;
275 
276  if(!a)
277  return ERROR_PTR;
278  if(n < 1 || m < 1)
279  return ERROR_MATRIX_SIZE;
280 
281  if(!a)
282  return ERROR_PTR;
283  if(n < 1 || m < 1)
284  return ERROR_SIZE;
285 
286  printf("\n%s = [ %% size [%d x %d] type: complex", name, n, m);
287 
288  for(p = 0; p < n; p++)
289  {
290  printf("\n");
291  for(q = 0; q < m; q++)
292  {
293  printf(format, RE(a[q*n + p]), IM(a[q*n + p]));
294  if(q == m-1)
295  printf(";");
296  else
297  printf(", ");
298  }
299  }
300  printf("];\n");
301 
302  return RES_OK;
303 }
304 
305 
306 
307 
308 
309 #ifdef DOXYGEN_ENGLISH
310 
311 #endif
312 #ifdef DOXYGEN_RUSSIAN
313 
314 #endif
315 int DSPL_API matrix_transpose(double* a, int n, int m, double* b)
316 {
317  int p, q, i, j, aind, bind;
318  if(!a || !b)
319  return ERROR_PTR;
320  if(n < 1 || m < 1)
321  return ERROR_MATRIX_SIZE;
322 
323 
324  for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK)
325  {
326  for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK)
327  {
328  for(i = 0; i < DSPL_MATRIX_BLOCK; i++)
329  {
330  for(j = 0; j < DSPL_MATRIX_BLOCK; j++)
331  {
332  aind = (q+j) * n + p + i;
333  bind = (p+i) * m + q + j;
334  b[bind] = a[aind];
335  }
336  }
337  }
338  }
339  for(i = p; i < n; i++)
340  for(j = 0; j < m; j++)
341  b[i*m + j] = a[j*n+i];
342 
343  for(i = 0; i < p; i++)
344  for(j = q; j < m; j++)
345  b[i*m + j] = a[j*n+i];
346 
347  return RES_OK;
348 }
349 
350 
351 
352 
353 
354 #ifdef DOXYGEN_ENGLISH
355 
356 #endif
357 #ifdef DOXYGEN_RUSSIAN
358 
359 #endif
360 int DSPL_API matrix_transpose_cmplx(complex_t* a, int n, int m, complex_t* b)
361 {
362  int p, q, i, j, aind, bind;
363 
364  if(!a || !b)
365  return ERROR_PTR;
366  if(n < 1 || m < 1)
367  return ERROR_MATRIX_SIZE;
368 
369  for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK)
370  {
371  for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK)
372  {
373  for(i = 0; i < DSPL_MATRIX_BLOCK; i++)
374  {
375  for(j = 0; j < DSPL_MATRIX_BLOCK; j++)
376  {
377  aind = (q+j) * n + p + i;
378  bind = (p+i) * m + q + j;
379  RE(b[bind]) = RE(a[aind]);
380  IM(b[bind]) = IM(a[aind]);
381  }
382  }
383  }
384  }
385  for(i = p; i < n; i++)
386  {
387  for(j = 0; j < m; j++)
388  {
389  RE(b[i*m + j]) = RE(a[j*n+i]);
390  IM(b[i*m + j]) = IM(a[j*n+i]);
391  }
392  }
393 
394  for(i = 0; i < p; i++)
395  {
396  for(j = q; j < m; j++)
397  {
398  RE(b[i*m + j]) = RE(a[j*n+i]);
399  IM(b[i*m + j]) = IM(a[j*n+i]);
400  }
401  }
402  return RES_OK;
403 }
404 
405 
406 
407 
408 
409 #ifdef DOXYGEN_ENGLISH
410 
411 #endif
412 #ifdef DOXYGEN_RUSSIAN
413 
414 #endif
415 int DSPL_API matrix_transpose_hermite(complex_t* a, int n, int m, complex_t* b)
416 {
417  int p, q, i, j, aind, bind;
418 
419  if(!a || !b)
420  return ERROR_PTR;
421  if(n < 1 || m < 1)
422  return ERROR_MATRIX_SIZE;
423 
424  for(p = 0; p < n - DSPL_MATRIX_BLOCK; p+=DSPL_MATRIX_BLOCK)
425  {
426  for(q = 0; q < m - DSPL_MATRIX_BLOCK; q+=DSPL_MATRIX_BLOCK)
427  {
428  for(i = 0; i < DSPL_MATRIX_BLOCK; i++)
429  {
430  for(j = 0; j < DSPL_MATRIX_BLOCK; j++)
431  {
432  aind = (q+j) * n + p + i;
433  bind = (p+i) * m + q + j;
434  RE(b[bind]) = RE(a[aind]);
435  IM(b[bind]) = -IM(a[aind]);
436  }
437  }
438  }
439  }
440  for(i = p; i < n; i++)
441  {
442  for(j = 0; j < m; j++)
443  {
444  RE(b[i*m + j]) = RE(a[j*n+i]);
445  IM(b[i*m + j]) = -IM(a[j*n+i]);
446  }
447  }
448 
449  for(i = 0; i < p; i++)
450  {
451  for(j = q; j < m; j++)
452  {
453  RE(b[i*m + j]) = RE(a[j*n+i]);
454  IM(b[i*m + j]) = -IM(a[j*n+i]);
455  }
456  }
457 
458  return RES_OK;
459 }
460 
461 
462 
463 
464 #ifdef DOXYGEN_ENGLISH
465 
466 #endif
467 #ifdef DOXYGEN_RUSSIAN
468 
469 #endif
470 int DSPL_API vector_dot(double* x, double* y, int n, double* p)
471 {
472  int inc = 1;
473 
474  if(!x || !y || !p)
475  return ERROR_PTR;
476  if(n<1)
477  return ERROR_SIZE;
478 
479  *p = ddot_(&n, x, &inc, y, &inc);
480 
481  return RES_OK;
482 }
483 
#define ERROR_LAPACK
Встроенная функция пакета LAPACK вернула код ошибки. Данная ошибка возвращается функцией,...
Definition: dspl.h:534
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:359
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:545
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:553
int matrix_eig_cmplx(complex_t *a, int n, complex_t *v, int *info)
Расчет собственных значений квадратной комплексной матрицы.
Definition: matrix.c:90
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:86
int matrix_eye(double *a, int n, int m)
Генерирование единичной вещественой матрицы размерности n x m.
Definition: matrix.c:153
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:497
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:417
#define ERROR_MATRIX_SIZE
Неверный размер матрицы.
Definition: dspl.h:537