GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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"
10static double tpi = 6.283185307179586;
11
12static void uortho(double *g, int n);
13
14double unfl();
15
16void 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
70static 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 r
float g
Definition: named_colr.c:8
Definition: ccmath.h:38
double re
Definition: ccmath.h:38
double im
Definition: ccmath.h:38
double unfl()
Definition: unfl.c:10
void unitary(Cpx *u, int n)
Definition: unitary.c:16