Loading...
Searching...
No Matches
unitary.c
Go to the documentation of this file.
/* unitary.c CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
#include <stdlib.h>
#include "ccmath.h"
static double tpi = 6.283185307179586;
static void uortho(double *g, int n);
double unfl(void);
void unitary(Cpx *u, int n)
{
int i, j, k, m;
Cpx h, *v, *e, *p, *r;
double *g, *q, a;
m = n * n;
g = (double *)calloc(n * n, sizeof(double));
v = (Cpx *)calloc(m + n, sizeof(Cpx));
e = v + m;
h.re = 1.;
h.im = 0.;
for (i = 0; i < n; ++i) {
a = tpi * unfl();
e[i].re = cos(a);
e[i].im = sin(a);
a = h.re * e[i].re - h.im * e[i].im;
h.im = h.im * e[i].re + h.re * e[i].im;
h.re = a;
}
h.im = -h.im;
for (i = 0; i < n; ++i) {
a = e[i].re * h.re - e[i].im * h.im;
e[i].im = e[i].re * h.im + e[i].im * h.re;
e[i].re = a;
}
uortho(g, n);
for (i = 0, p = v, q = g; i < n; ++i) {
for (j = 0; j < n; ++j)
(p++)->re = *q++;
}
for (i = 0, p = v; i < n; ++i) {
for (j = 0, h = e[i]; j < n; ++j, ++p) {
a = h.re * p->re - h.im * p->im;
p->im = h.im * p->re + h.re * p->im;
p->re = a;
}
}
uortho(g, n);
for (i = m = 0, p = u; i < n; ++i, m += n) {
for (j = 0; j < n; ++j, ++p) {
p->re = p->im = 0.;
for (k = 0, q = g + m, r = v + j; k < n; ++k, r += n) {
p->re += *q * r->re;
p->im += *q++ * r->im;
}
}
}
free(g);
free(v);
}
static void uortho(double *g, int n)
{
int i, j, k, m;
double *p, *q, c, s, a;
for (i = 0, p = g; i < n; ++i) {
for (j = 0; j < n; ++j) {
if (i == j)
*p++ = 1.;
else
*p++ = 0.;
}
}
for (i = 0, m = n - 1; i < m; ++i) {
for (j = i + 1; j < n; ++j) {
a = tpi * unfl();
c = cos(a);
s = sin(a);
p = g + n * i;
q = g + n * j;
for (k = 0; k < n; ++k) {
a = *p * c + *q * s;
*q = *q * c - *p * s;
*p++ = a;
++q;
}
}
}
}
ccmath.h
Cpx
struct complex Cpx
Definition ccmath.h:41
r
double r
Definition driver/set_window.c:5
g
float g
Definition named_colr.c:7
complex::re
double re
Definition ccmath.h:39
complex::im
double im
Definition ccmath.h:39
unfl
double unfl(void)
Definition unfl.c:10
unitary
void unitary(Cpx *u, int n)
Definition unitary.c:18