1 /* matrix.c: Matrix operations
2
3 Copyright 2003, 2005 Bjoern Butscher, Hendrik Weimer
4
5 This file is part of libquantum
6
7 libquantum is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published
9 by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 libquantum is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with libquantum; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
20 MA 02110-1301, USA
21
22 */
23
24 #include <stdlib.h>
25 #include <stdio.h>
26
27 #include "matrix.h"
28 #include "config.h"
29 #include "complex.h"
30 #include "error.h"
31
32 /* Statistics of the memory consumption */
33
quantum_memman(long change)34 unsigned long quantum_memman(long change)
35 {
36 static long mem = 0, max = 0;
37
38 mem += change;
39
40 if(mem > max)
41 max = mem;
42
43 return mem;
44 }
45
46 /* Create a new COLS x ROWS matrix */
47
48 quantum_matrix
quantum_new_matrix(int cols,int rows)49 quantum_new_matrix(int cols, int rows)
50 {
51 quantum_matrix m;
52
53 m.rows = rows;
54 m.cols = cols;
55 m.t = calloc(cols * rows, sizeof(COMPLEX_FLOAT));
56
57 #if (DEBUG_MEM)
58 printf("allocating %i bytes of memory for %ix%i matrix at 0x%X\n",
59 sizeof(COMPLEX_FLOAT) * cols * rows, cols, rows, (int) m.t);
60 #endif
61
62 if(!m.t)
63 quantum_error(QUANTUM_ENOMEM);
64
65 quantum_memman(sizeof(COMPLEX_FLOAT) * cols * rows);
66
67 return m;
68 }
69
70 /* Delete a matrix */
71
72 void
quantum_delete_matrix(quantum_matrix * m)73 quantum_delete_matrix(quantum_matrix *m)
74 {
75 #if (DEBUG_MEM)
76 printf("freeing %i bytes of memory for %ix%i matrix at 0x%X\n",
77 sizeof(COMPLEX_FLOAT) * m->cols * m->rows, m->cols, m->rows,
78 (int) m->t);
79 #endif
80
81 free(m->t);
82 quantum_memman(-sizeof(COMPLEX_FLOAT) * m->cols * m->rows);
83 m->t=0;
84 }
85
86 /* Print the contents of a matrix to stdout */
87
88 void
quantum_print_matrix(quantum_matrix m)89 quantum_print_matrix(quantum_matrix m)
90 {
91 int i, j, z=0;
92 int print_imag = 0;
93 /* int l; */
94
95 for(i=0; i<m.rows; i++)
96 {
97 for(j=0; j<m.cols; j++)
98 {
99 if(quantum_imag(M(m, j, i))/quantum_real(M(m, j, i)) > 1e-3)
100 print_imag = 1;
101 }
102 }
103
104 while ((1 << z++) < m.rows);
105 z--;
106
107 for(i=0; i<m.rows; i++)
108 {
109 /* for (l=z-1; l>=0; l--)
110 {
111 if ((l % 4 == 3))
112 printf(" ");
113 printf("%i", (i >> l) & 1);
114 } */
115
116 for(j=0; j<m.cols; j++)
117 {
118 if(print_imag)
119 printf("%3.3f%+.3fi ", quantum_real(M(m, j, i)),
120 quantum_imag(M(m, j, i)));
121 else
122 // printf("%3.3f ", quantum_real(M(m, j, i)));
123 printf("%+.1f ", quantum_real(M(m, j, i)));
124 }
125
126 printf("\n");
127 }
128 printf("\n");
129 }
130
131 /* Matrix multiplication */
132
quantum_mmult(quantum_matrix A,quantum_matrix B)133 quantum_matrix quantum_mmult(quantum_matrix A, quantum_matrix B)
134 {
135 int i, j, k;
136 quantum_matrix C;
137
138 if(A.cols != B.rows)
139 quantum_error(QUANTUM_EMSIZE);
140
141 C = quantum_new_matrix(B.cols, A.rows);
142
143 for(i=0; i<B.cols; i++)
144 {
145 for(j=0; j<A.rows; j++)
146 {
147 for(k=0; k<B.rows; k++)
148 M(C, i, j) += M(A, k, j) * M(B, i, k);
149 }
150 }
151
152 return C;
153 }
154
155 /* Compute the adjoint of a matrix */
156
157 void
quantum_adjoint(quantum_matrix * m)158 quantum_adjoint(quantum_matrix *m)
159 {
160 int i, j;
161 COMPLEX_FLOAT tmp;
162 quantum_matrix A = *m;
163
164 for(i=0; i<m->cols; i++)
165 {
166 for(j=0;j<i;j++)
167 {
168 tmp = M(A, i, j);
169 M(A, i, j) = quantum_conj(M(A, j, i));
170 M(A, j, i) = quantum_conj(tmp);
171 }
172 }
173 }
174
175