GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
9#include <stdlib.h>
10#include "ccmath.h"
11
12static double tpi = 6.283185307179586;
13
14static void uortho(double *g, int n);
15
16double unfl(void);
17
18void unitary(Cpx *u, int n)
19{
20 int i, j, k, m;
21
22 Cpx h, *v, *e, *p, *r;
23
24 double *g, *q, a;
25
26 m = n * n;
27 g = (double *)calloc(n * n, sizeof(double));
28 v = (Cpx *)calloc(m + n, sizeof(Cpx));
29 e = v + m;
30 h.re = 1.;
31 h.im = 0.;
32 for (i = 0; i < n; ++i) {
33 a = tpi * unfl();
34 e[i].re = cos(a);
35 e[i].im = sin(a);
36 a = h.re * e[i].re - h.im * e[i].im;
37 h.im = h.im * e[i].re + h.re * e[i].im;
38 h.re = a;
39 }
40 h.im = -h.im;
41 for (i = 0; i < n; ++i) {
42 a = e[i].re * h.re - e[i].im * h.im;
43 e[i].im = e[i].re * h.im + e[i].im * h.re;
44 e[i].re = a;
45 }
46 uortho(g, n);
47 for (i = 0, p = v, q = g; i < n; ++i) {
48 for (j = 0; j < n; ++j)
49 (p++)->re = *q++;
50 }
51 for (i = 0, p = v; i < n; ++i) {
52 for (j = 0, h = e[i]; j < n; ++j, ++p) {
53 a = h.re * p->re - h.im * p->im;
54 p->im = h.im * p->re + h.re * p->im;
55 p->re = a;
56 }
57 }
58 uortho(g, n);
59 for (i = m = 0, p = u; i < n; ++i, m += n) {
60 for (j = 0; j < n; ++j, ++p) {
61 p->re = p->im = 0.;
62 for (k = 0, q = g + m, r = v + j; k < n; ++k, r += n) {
63 p->re += *q * r->re;
64 p->im += *q++ * r->im;
65 }
66 }
67 }
68 free(g);
69 free(v);
70}
71
72static void uortho(double *g, int n)
73{
74 int i, j, k, m;
75
76 double *p, *q, c, s, a;
77
78 for (i = 0, p = g; i < n; ++i) {
79 for (j = 0; j < n; ++j) {
80 if (i == j)
81 *p++ = 1.;
82 else
83 *p++ = 0.;
84 }
85 }
86 for (i = 0, m = n - 1; i < m; ++i) {
87 for (j = i + 1; j < n; ++j) {
88 a = tpi * unfl();
89 c = cos(a);
90 s = sin(a);
91 p = g + n * i;
92 q = g + n * j;
93 for (k = 0; k < n; ++k) {
94 a = *p * c + *q * s;
95 *q = *q * c - *p * s;
96 *p++ = a;
97 ++q;
98 }
99 }
100 }
101}
struct complex Cpx
Definition ccmath.h:41
double r
float g
Definition named_colr.c:7
double re
Definition ccmath.h:39
double im
Definition ccmath.h:39
double unfl(void)
Definition unfl.c:10
void unitary(Cpx *u, int n)
Definition unitary.c:18