GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
unitary.c
Go to the documentation of this file.
1 /* unitary.c CCMATH mathematics library source code.
2  *
3  * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
4  * This code may be redistributed under the terms of the GNU library
5  * public license (LGPL). ( See the lgpl.license file for details.)
6  * ------------------------------------------------------------------------
7  */
8 #include <stdlib.h>
9 #include "ccmath.h"
10 static double tpi = 6.283185307179586;
11 
12 static void uortho(double *g, int n);
13 
14 double unfl();
15 
16 void unitary(Cpx * u, int n)
17 {
18  int i, j, k, m;
19 
20  Cpx h, *v, *e, *p, *r;
21 
22  double *g, *q, a;
23 
24  m = n * n;
25  g = (double *)calloc(n * n, sizeof(double));
26  v = (Cpx *) calloc(m + n, sizeof(Cpx));
27  e = v + m;
28  h.re = 1.;
29  h.im = 0.;
30  for (i = 0; i < n; ++i) {
31  a = tpi * unfl();
32  e[i].re = cos(a);
33  e[i].im = sin(a);
34  a = h.re * e[i].re - h.im * e[i].im;
35  h.im = h.im * e[i].re + h.re * e[i].im;
36  h.re = a;
37  }
38  h.im = -h.im;
39  for (i = 0; i < n; ++i) {
40  a = e[i].re * h.re - e[i].im * h.im;
41  e[i].im = e[i].re * h.im + e[i].im * h.re;
42  e[i].re = a;
43  }
44  uortho(g, n);
45  for (i = 0, p = v, q = g; i < n; ++i) {
46  for (j = 0; j < n; ++j)
47  (p++)->re = *q++;
48  }
49  for (i = 0, p = v; i < n; ++i) {
50  for (j = 0, h = e[i]; j < n; ++j, ++p) {
51  a = h.re * p->re - h.im * p->im;
52  p->im = h.im * p->re + h.re * p->im;
53  p->re = a;
54  }
55  }
56  uortho(g, n);
57  for (i = m = 0, p = u; i < n; ++i, m += n) {
58  for (j = 0; j < n; ++j, ++p) {
59  p->re = p->im = 0.;
60  for (k = 0, q = g + m, r = v + j; k < n; ++k, r += n) {
61  p->re += *q * r->re;
62  p->im += *q++ * r->im;
63  }
64  }
65  }
66  free(g);
67  free(v);
68 }
69 
70 static void uortho(double *g, int n)
71 {
72  int i, j, k, m;
73 
74  double *p, *q, c, s, a;
75 
76  for (i = 0, p = g; i < n; ++i) {
77  for (j = 0; j < n; ++j) {
78  if (i == j)
79  *p++ = 1.;
80  else
81  *p++ = 0.;
82  }
83  }
84  for (i = 0, m = n - 1; i < m; ++i) {
85  for (j = i + 1; j < n; ++j) {
86  a = tpi * unfl();
87  c = cos(a);
88  s = sin(a);
89  p = g + n * i;
90  q = g + n * j;
91  for (k = 0; k < n; ++k) {
92  a = *p * c + *q * s;
93  *q = *q * c - *p * s;
94  *p++ = a;
95  ++q;
96  }
97  }
98  }
99 }
double unfl()
Definition: unfl.c:10
void unitary(Cpx *u, int n)
Definition: unitary.c:16
void free(void *)
float g
Definition: named_colr.c:8
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38
Definition: la.h:54
double r
Definition: r_raster.c:39