8 #include "exponentialsum.h"
9 #include "shrinkstar.h"
11 #include "measurepauli.h"
12 #include "innerproduct.h"
13 #include "randomstabilizerstate.h"
15 #define ZEROTHRESHOLD (0.00000001)
17 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
19 // order of matrix elements is [row][column]!!!
21 int main(int argc, char *argv[])
25 printf("strongsim12_rellerr argument: \"number of stabilizer state samples\"\n");
29 int NUMSTABSTATESAMPLES = atoi(argv[1]); // number of stabilizer state samples
31 int N; // number of qubits
35 printf("'N' needs to be a multiple of 12 for a k=12 tensor factor decomposition!\n");
39 int T; // number of T gate magic states (set to the first 'K' of the 'N' qubits -- the rest are set to the '0' computational basis state)
42 int omega[N]; // max of N measurements
43 int alpha[N][N], beta[N][N], gamma[N][N], delta[N][N]; // max of N measurements of N Paulis
51 srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
53 fp = fopen("Pd.txt", "r");
55 if(fscanf(fp, "%s", buff) == EOF) {
56 printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
64 Pd = calloc(PdN, sizeof(double*));
66 Pd[i] = calloc(PdN+1, sizeof(double));
70 for(i=1; i<PdN; i++) {
73 if(fscanf(fp, "%s", buff) == EOF) {
74 printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
77 Pd[i][j] = atof(buff);
78 //printf("%e ", Pd[i][j]);
82 //printf("total=%f\n", tmp);
86 double complex coeffb60 = (-16.0+12.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0);
87 double complex coeffb66 = (96.0-68.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)/8.0*cpow(2.0,3.0);
88 double complex coeffe6 = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(6.0*PI*I/8.0)*cpow(2.0,2.5);
89 double complex coeffo6 = (-14.0+10.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-14.0*PI*I/8.0)*cpow(2.0,2.5);
90 double complex coeffk6 = (7.0-5.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(-8.0*PI*I/8.0)*4.0*csqrt(2.0)*cpow(2.0,0.5);
91 double complex coeffphiprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5);
92 double complex coeffphidprime = (10.0-7.0*sqrt(2.0)+0.0*I)*pow(cos(PI/8.0),6)*cexp(2.0*PI*I/8.0)*cpow(2.0,2.5);
95 int nn1 = 6; int kk1 = 6; int (*(GG1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GGBar1[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int hh1[] = {0, 0, 0, 0, 0, 0}; int QQ1 = 0; int DD1[] = {0, 0, 0, 0, 0, 0}; int (*(JJ1[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} };
97 int nn2 = 6; int kk2 = 6; int (*(GG2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int (*(GGBar2[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}}; int hh2[] = {0, 0, 0, 0, 0, 0}; int QQ2 = 4; int DD2[] = {4, 4, 4, 4, 4, 4}; int (*(JJ2[])) = { (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0} };
99 int nn3 = 6; int kk3 = 5; int (*(GG3[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar3[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh3[] = {1, 0, 0, 0, 0, 0}; int QQ3 = 4; int DD3[] = {0, 0, 0, 0, 0}; int (*(JJ3[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} };
101 int nn4 = 6; int kk4 = 5; int (*(GG4[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar4[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh4[] = {0, 0, 0, 0, 0, 0}; int QQ4 = 4; int DD4[] = {4, 4, 4, 4, 4}; int (*(JJ4[])) = { (int[]) {0, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 0} };
103 int nn5 = 6; int kk5 = 1; int (*(GG5[])) = { (int[]) {1, 1, 1, 1, 1, 1}, (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1} }; int (*(GGBar5[])) = { (int[]) {1, 0, 0, 0, 0, 0}, (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1} }; int hh5[] = {1, 1, 1, 1, 1, 1}; int QQ5 = 6; int DD5[] = {2}; int (*(JJ5[])) = { (int[]) {4} };
105 int nn6 = 6; int kk6 = 5; int (*(GG6[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar6[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh6[] = {1, 0, 0, 0, 0, 0}; int QQ6 = 0; int DD6[] = {0, 0, 0, 0, 0}; int (*(JJ6[])) = { (int[]) {0, 4, 0, 0, 4}, (int[]) {4, 0, 4, 0, 0}, (int[]) {0, 4, 0, 4, 0}, (int[]) {0, 0, 4, 0, 4}, (int[]) {4, 0, 0, 4, 0} };
107 int nn7 = 6; int kk7 = 5; int (*(GG7[])) = { (int[]) {1, 1, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0} }; int (*(GGBar7[])) = { (int[]) {0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1} }; int hh7[] = {1, 0, 0, 0, 0, 0}; int QQ7 = 0; int DD7[] = {0, 0, 0, 0, 0}; int (*(JJ7[])) = { (int[]) {0, 0, 4, 4, 0}, (int[]) {0, 0, 0, 4, 4}, (int[]) {4, 0, 0, 0, 4}, (int[]) {4, 4, 0, 0, 0}, (int[]) {0, 4, 4, 0, 0} };
110 int n1 = nn1+nn1; int k1 = kk1+kk1; int **G1; int **GBar1; int *h1; int Q1; int *D1; int **J1;
111 G1 = calloc(nn1+nn1,sizeof(int*));
112 GBar1 = calloc(nn1+nn1,sizeof(int*));
113 h1 = calloc(nn1+nn1,sizeof(int));
114 D1 = calloc(kk1+kk1,sizeof(int));
115 J1 = calloc(kk1+kk1,sizeof(int*));
116 for(i=0; i<nn1+nn1; i++) {
117 G1[i] = calloc(nn1+nn1, sizeof(int)); GBar1[i] = calloc(nn1+nn1, sizeof(int));
119 for(i=0; i<kk1+kk1; i++)
120 J1[i] = calloc(kk1+kk1, sizeof(int));
121 //appendBlockMatrix(GG1, GG1, G1, nn1, nn1, nn1, nn1);
122 addSubMatrix(GG1, G1, kk1, nn1, 0, 0, 0, 0);
123 addSubMatrix(GG1, G1, kk1, nn1, 0, 0, kk1, nn1);
124 addSubMatrix(GG1, G1, nn1-kk1, nn1, kk1, 0, kk1+kk1, 0); // put the out-of-subspace vectors at the end
125 addSubMatrix(GG1, G1, nn1-kk1, nn1, kk1, 0, kk1+kk1+nn1-kk1, nn1);
126 //appendBlockMatrix(GGBar1, GGBar1, GBar1, nn1, nn1, nn1, nn1);
127 addSubMatrix(GGBar1, GBar1, kk1, nn1, 0, 0, 0, 0);
128 addSubMatrix(GGBar1, GBar1, kk1, nn1, 0, 0, kk1, nn1);
129 addSubMatrix(GGBar1, GBar1, nn1-kk1, nn1, kk1, 0, kk1+kk1, 0); // put the out-of-subspace vectors at the end
130 addSubMatrix(GGBar1, GBar1, nn1-kk1, nn1, kk1, 0, kk1+kk1+nn1-kk1, nn1);
131 appendVector(hh1, hh1, h1, nn1, nn1);
133 appendVector(DD1, DD1, D1, kk1, kk1);
134 appendBlockMatrix(JJ1, JJ1, J1, kk1, kk1, kk1, kk1);
135 //addSubMatrix(JJ1, J1, kk1, kk1, 0, 0, 0, 0);
136 //addSubMatrix(JJ1, J1, kk1, kk1, 0, 0, kk1, kk1);
139 /*int n2 = nn1+nn2; int k2 = kk1+kk2; int **G2; int **GBar2; int *h2; int Q2; int *D2; int **J2;
140 G2 = calloc(nn1+nn2,sizeof(int*));
141 GBar2 = calloc(nn1+nn2,sizeof(int*));
142 h2 = calloc(nn1+nn2,sizeof(int));
143 D2 = calloc(kk1+kk2,sizeof(int));
144 J2 = calloc(kk1+kk2,sizeof(int*));
145 for(i=0; i<nn1+nn2; i++) {
146 G2[i] = calloc(nn1+nn2, sizeof(int)); GBar2[i] = calloc(nn1+nn2, sizeof(int));
148 for(i=0; i<kk1+kk2; i++)
149 J2[i] = calloc(kk1+kk2, sizeof(int));
150 //appendBlockMatrix(GG2, GG1, G2, nn2, nn2, nn1, nn1);
151 addSubMatrix(GG2, G2, kk2, nn2, 0, 0, 0, 0);
152 addSubMatrix(GG1, G2, kk1, nn1, 0, 0, kk2, nn2);
153 addSubMatrix(GG2, G2, nn2-kk2, nn2, kk2, 0, kk2+kk1, 0); // put the out-of-subspace vectors at the end
154 addSubMatrix(GG1, G2, nn1-kk1, nn1, kk1, 0, kk2+kk1+nn2-kk2, nn1);
155 //appendBlockMatrix(GGBar2, GGBar1, GBar2, nn2, nn2, nn1, nn1);
156 addSubMatrix(GGBar2, GBar2, kk2, nn2, 0, 0, 0, 0);
157 addSubMatrix(GGBar1, GBar2, kk1, nn1, 0, 0, kk2, nn2);
158 addSubMatrix(GGBar2, GBar2, nn2-kk2, nn2, kk2, 0, kk2+kk1, 0); // put the out-of-subspace vectors at the end
159 addSubMatrix(GGBar1, GBar2, nn1-kk1, nn1, kk1, 0, kk2+kk1+nn2-kk2, nn1);
160 appendVector(hh2, hh1, h2, nn2, nn1);
162 appendVector(DD2, DD1, D2, kk2, kk1);
163 appendBlockMatrix(JJ2, JJ1, J2, kk2, kk2, kk1, kk1);*/
166 int n2 = nn1+nn3; int k2 = kk1+kk3; int **G2; int **GBar2; int *h2; int Q2; int *D2; int **J2;
167 G2 = calloc(nn1+nn3,sizeof(int*));
168 GBar2 = calloc(nn1+nn3,sizeof(int*));
169 h2 = calloc(nn1+nn3,sizeof(int));
170 D2 = calloc(kk1+kk3,sizeof(int));
171 J2 = calloc(kk1+kk3,sizeof(int*));
172 for(i=0; i<nn1+nn3; i++) {
173 G2[i] = calloc(nn1+nn3, sizeof(int)); GBar2[i] = calloc(nn1+nn3, sizeof(int));
175 for(i=0; i<kk1+kk3; i++)
176 J2[i] = calloc(kk1+kk3, sizeof(int));
177 //appendBlockMatrix(GG3, GG1, G3, nn3, nn3, nn1, nn1);
178 addSubMatrix(GG3, G2, kk3, nn3, 0, 0, 0, 0);
179 addSubMatrix(GG1, G2, kk1, nn1, 0, 0, kk3, nn3);
180 addSubMatrix(GG3, G2, nn3-kk3, nn3, kk3, 0, kk3+kk1, 0); // put the out-of-subspace vectors at the end
181 addSubMatrix(GG1, G2, nn1-kk1, nn1, kk1, 0, kk3+kk1+nn3-kk3, nn1);
182 //appendBlockMatrix(GGBar3, GGBar1, GBar3, nn3, nn3, nn1, nn1);
183 addSubMatrix(GGBar3, GBar2, kk3, nn3, 0, 0, 0, 0);
184 addSubMatrix(GGBar1, GBar2, kk1, nn1, 0, 0, kk3, nn3);
185 addSubMatrix(GGBar3, GBar2, nn3-kk3, nn3, kk3, 0, kk3+kk1, 0); // put the out-of-subspace vectors at the end
186 addSubMatrix(GGBar1, GBar2, nn1-kk1, nn1, kk1, 0, kk3+kk1+nn3-kk3, nn1);
187 appendVector(hh3, hh1, h2, nn3, nn1);
189 appendVector(DD3, DD1, D2, kk3, kk1);
190 appendBlockMatrix(JJ3, JJ1, J2, kk3, kk3, kk1, kk1);
193 int n3 = nn1+nn4; int k3 = kk1+kk4; int **G3; int **GBar3; int *h3; int Q3; int *D3; int **J3;
194 G3 = calloc(nn1+nn4,sizeof(int*));
195 GBar3 = calloc(nn1+nn4,sizeof(int*));
196 h3 = calloc(nn1+nn4,sizeof(int));
197 D3 = calloc(kk1+kk4,sizeof(int));
198 J3 = calloc(kk1+kk4,sizeof(int*));
199 for(i=0; i<nn1+nn4; i++) {
200 G3[i] = calloc(nn1+nn4, sizeof(int)); GBar3[i] = calloc(nn1+nn4, sizeof(int));
202 for(i=0; i<kk1+kk4; i++)
203 J3[i] = calloc(kk1+kk4, sizeof(int));
204 //appendBlockMatrix(GG4, GG1, G4, nn4, nn4, nn1, nn1);
205 addSubMatrix(GG4, G3, kk4, nn4, 0, 0, 0, 0);
206 addSubMatrix(GG1, G3, kk1, nn1, 0, 0, kk4, nn4);
207 addSubMatrix(GG4, G3, nn4-kk4, nn4, kk4, 0, kk4+kk1, 0); // put the out-of-subspace vectors at the end
208 addSubMatrix(GG1, G3, nn1-kk1, nn1, kk1, 0, kk4+kk1+nn4-kk4, nn1);
209 //appendBlockMatrix(GGBar4, GGBar1, GBar4, nn4, nn4, nn1, nn1);
210 addSubMatrix(GGBar4, GBar3, kk4, nn4, 0, 0, 0, 0);
211 addSubMatrix(GGBar1, GBar3, kk1, nn1, 0, 0, kk4, nn4);
212 addSubMatrix(GGBar4, GBar3, nn4-kk4, nn4, kk4, 0, kk4+kk1, 0); // put the out-of-subspace vectors at the end
213 addSubMatrix(GGBar1, GBar3, nn1-kk1, nn1, kk1, 0, kk4+kk1+nn4-kk4, nn1);
214 appendVector(hh4, hh1, h3, nn4, nn1);
216 appendVector(DD4, DD1, D3, kk4, kk1);
217 appendBlockMatrix(JJ4, JJ1, J3, kk4, kk4, kk1, kk1);
220 int n4 = nn1+nn5; int k4 = kk1+kk5; int **G4; int **GBar4; int *h4; int Q4; int *D4; int **J4;
221 G4 = calloc(nn1+nn5,sizeof(int*));
222 GBar4 = calloc(nn1+nn5,sizeof(int*));
223 h4 = calloc(nn1+nn5,sizeof(int));
224 D4 = calloc(kk1+kk5,sizeof(int));
225 J4 = calloc(kk1+kk5,sizeof(int*));
226 for(i=0; i<nn1+nn5; i++) {
227 G4[i] = calloc(nn1+nn5, sizeof(int)); GBar4[i] = calloc(nn1+nn5, sizeof(int));
229 for(i=0; i<kk1+kk5; i++)
230 J4[i] = calloc(kk1+kk5, sizeof(int));
231 //appendBlockMatrix(GG5, GG1, G5, nn5, nn5, nn1, nn1);
232 addSubMatrix(GG5, G4, kk5, nn5, 0, 0, 0, 0);
233 addSubMatrix(GG1, G4, kk1, nn1, 0, 0, kk5, nn5);
234 addSubMatrix(GG5, G4, nn5-kk5, nn5, kk5, 0, kk5+kk1, 0); // put the out-of-subspace vectors at the end
235 addSubMatrix(GG1, G4, nn1-kk1, nn1, kk1, 0, kk5+kk1+nn5-kk5, nn1);
236 //appendBlockMatrix(GGBar5, GGBar1, GBar5, nn5, nn5, nn1, nn1);
237 addSubMatrix(GGBar5, GBar4, kk5, nn5, 0, 0, 0, 0);
238 addSubMatrix(GGBar1, GBar4, kk1, nn1, 0, 0, kk5, nn5);
239 addSubMatrix(GGBar5, GBar4, nn5-kk5, nn5, kk5, 0, kk5+kk1, 0); // put the out-of-subspace vectors at the end
240 addSubMatrix(GGBar1, GBar4, nn1-kk1, nn1, kk1, 0, kk5+kk1+nn5-kk5, nn1);
241 appendVector(hh5, hh1, h4, nn5, nn1);
243 appendVector(DD5, DD1, D4, kk5, kk1);
244 appendBlockMatrix(JJ5, JJ1, J4, kk5, kk5, kk1, kk1);
247 int n5 = nn1+nn6; int k5 = kk1+kk6; int **G5; int **GBar5; int *h5; int Q5; int *D5; int **J5;
248 G5 = calloc(nn1+nn6,sizeof(int*));
249 GBar5 = calloc(nn1+nn6,sizeof(int*));
250 h5 = calloc(nn1+nn6,sizeof(int));
251 D5 = calloc(kk1+kk6,sizeof(int));
252 J5 = calloc(kk1+kk6,sizeof(int*));
253 for(i=0; i<nn1+nn6; i++) {
254 G5[i] = calloc(nn1+nn6, sizeof(int)); GBar5[i] = calloc(nn1+nn6, sizeof(int));
256 for(i=0; i<kk1+kk6; i++)
257 J5[i] = calloc(kk1+kk6, sizeof(int));
258 //appendBlockMatrix(GG6, GG1, G6, nn6, nn6, nn1, nn1);
259 addSubMatrix(GG6, G5, kk6, nn6, 0, 0, 0, 0);
260 addSubMatrix(GG1, G5, kk1, nn1, 0, 0, kk6, nn6);
261 addSubMatrix(GG6, G5, nn6-kk6, nn6, kk6, 0, kk6+kk1, 0); // put the out-of-subspace vectors at the end
262 addSubMatrix(GG1, G5, nn1-kk1, nn1, kk1, 0, kk6+kk1+nn6-kk6, nn1);
263 //appendBlockMatrix(GGBar6, GGBar1, GBar6, nn6, nn6, nn1, nn1);
264 addSubMatrix(GGBar6, GBar5, kk6, nn6, 0, 0, 0, 0);
265 addSubMatrix(GGBar1, GBar5, kk1, nn1, 0, 0, kk6, nn6);
266 addSubMatrix(GGBar6, GBar5, nn6-kk6, nn6, kk6, 0, kk6+kk1, 0); // put the out-of-subspace vectors at the end
267 addSubMatrix(GGBar1, GBar5, nn1-kk1, nn1, kk1, 0, kk6+kk1+nn6-kk6, nn1);
268 appendVector(hh6, hh1, h5, nn6, nn1);
270 appendVector(DD6, DD1, D5, kk6, kk1);
271 appendBlockMatrix(JJ6, JJ1, J5, kk6, kk6, kk1, kk1);
274 int n6 = nn1+nn7; int k6 = kk1+kk7; int **G6; int **GBar6; int *h6; int Q6; int *D6; int **J6;
275 G6 = calloc(nn1+nn7,sizeof(int*));
276 GBar6 = calloc(nn1+nn7,sizeof(int*));
277 h6 = calloc(nn1+nn7,sizeof(int));
278 D6 = calloc(kk1+kk7,sizeof(int));
279 J6 = calloc(kk1+kk7,sizeof(int*));
280 for(i=0; i<nn1+nn7; i++) {
281 G6[i] = calloc(nn1+nn7, sizeof(int)); GBar6[i] = calloc(nn1+nn7, sizeof(int));
283 for(i=0; i<kk1+kk7; i++)
284 J6[i] = calloc(kk1+kk7, sizeof(int));
285 //appendBlockMatrix(GG7, GG1, G7, nn7, nn7, nn1, nn1);
286 addSubMatrix(GG7, G6, kk7, nn7, 0, 0, 0, 0);
287 addSubMatrix(GG1, G6, kk1, nn1, 0, 0, kk7, nn7);
288 addSubMatrix(GG7, G6, nn7-kk7, nn7, kk7, 0, kk7+kk1, 0); // put the out-of-subspace vectors at the end
289 addSubMatrix(GG1, G6, nn1-kk1, nn1, kk1, 0, kk7+kk1+nn7-kk7, nn1);
290 //appendBlockMatrix(GGBar7, GGBar1, GBar7, nn7, nn7, nn1, nn1);
291 addSubMatrix(GGBar7, GBar6, kk7, nn7, 0, 0, 0, 0);
292 addSubMatrix(GGBar1, GBar6, kk1, nn1, 0, 0, kk7, nn7);
293 addSubMatrix(GGBar7, GBar6, nn7-kk7, nn7, kk7, 0, kk7+kk1, 0); // put the out-of-subspace vectors at the end
294 addSubMatrix(GGBar1, GBar6, nn1-kk1, nn1, kk1, 0, kk7+kk1+nn7-kk7, nn1);
295 appendVector(hh7, hh1, h6, nn7, nn1);
297 appendVector(DD7, DD1, D6, kk7, kk1);
298 appendBlockMatrix(JJ7, JJ1, J6, kk7, kk7, kk1, kk1);
301 /*int n8 = nn2+nn1; int k8 = kk2+kk1; int **G8; int **GBar8; int *h8; int Q8; int *D8; int **J8;
302 G8 = calloc(nn2+nn1,sizeof(int*));
303 GBar8 = calloc(nn2+nn1,sizeof(int*));
304 h8 = calloc(nn2+nn1,sizeof(int));
305 D8 = calloc(kk2+kk1,sizeof(int));
306 J8 = calloc(kk2+kk1,sizeof(int*));
307 for(i=0; i<nn2+nn1; i++) {
308 G8[i] = calloc(nn2+nn1, sizeof(int)); GBar8[i] = calloc(nn2+nn1, sizeof(int));
310 for(i=0; i<kk2+kk1; i++)
311 J8[i] = calloc(kk2+kk1, sizeof(int));
312 //appendBlockMatrix(GG1, GG2, G8, nn1, nn1, nn2, nn2);
313 addSubMatrix(GG1, G8, kk1, nn1, 0, 0, 0, 0);
314 addSubMatrix(GG2, G8, kk2, nn2, 0, 0, kk1, nn1);
315 addSubMatrix(GG1, G8, nn1-kk1, nn1, kk1, 0, kk1+kk2, 0); // put the out-of-subspace vectors at the end
316 addSubMatrix(GG2, G8, nn2-kk2, nn2, kk2, 0, kk1+kk2+nn1-kk1, nn2);
317 //appendBlockMatrix(GGBar1, GGBar2, GBar8, nn1, nn1, nn2, nn2);
318 addSubMatrix(GGBar1, GBar8, kk1, nn1, 0, 0, 0, 0);
319 addSubMatrix(GGBar2, GBar8, kk2, nn2, 0, 0, kk1, nn1);
320 addSubMatrix(GGBar1, GBar8, nn1-kk1, nn1, kk1, 0, kk1+kk2, 0); // put the out-of-subspace vectors at the end
321 addSubMatrix(GGBar2, GBar8, nn2-kk2, nn2, kk2, 0, kk1+kk2+nn1-kk1, nn2);
322 appendVector(hh1, hh2, h8, nn1, nn2);
324 appendVector(DD1, DD2, D8, kk1, kk2);
325 appendBlockMatrix(JJ1, JJ2, J8, kk1, kk1, kk2, kk2);*/
328 int n7 = nn2+nn2; int k7 = kk2+kk2; int **G7; int **GBar7; int *h7; int Q7; int *D7; int **J7;
329 G7 = calloc(nn2+nn2,sizeof(int*));
330 GBar7 = calloc(nn2+nn2,sizeof(int*));
331 h7 = calloc(nn2+nn2,sizeof(int));
332 D7 = calloc(kk2+kk2,sizeof(int));
333 J7 = calloc(kk2+kk2,sizeof(int*));
334 for(i=0; i<nn2+nn2; i++) {
335 G7[i] = calloc(nn2+nn2, sizeof(int)); GBar7[i] = calloc(nn2+nn2, sizeof(int));
337 for(i=0; i<kk2+kk2; i++)
338 J7[i] = calloc(kk2+kk2, sizeof(int));
339 //appendBlockMatrix(GG2, GG2, G9, nn2, nn2, nn2, nn2);
340 addSubMatrix(GG2, G7, kk2, nn2, 0, 0, 0, 0);
341 addSubMatrix(GG2, G7, kk2, nn2, 0, 0, kk2, nn2);
342 addSubMatrix(GG2, G7, nn2-kk2, nn2, kk2, 0, kk2+kk2, 0); // put the out-of-subspace vectors at the end
343 addSubMatrix(GG2, G7, nn2-kk2, nn2, kk2, 0, kk2+kk2+nn2-kk2, nn2);
344 //appendBlockMatrix(GGBar2, GGBar2, GBar9, nn2, nn2, nn2, nn2);
345 addSubMatrix(GGBar2, GBar7, kk2, nn2, 0, 0, 0, 0);
346 addSubMatrix(GGBar2, GBar7, kk2, nn2, 0, 0, kk2, nn2);
347 addSubMatrix(GGBar2, GBar7, nn2-kk2, nn2, kk2, 0, kk2+kk2, 0); // put the out-of-subspace vectors at the end
348 addSubMatrix(GGBar2, GBar7, nn2-kk2, nn2, kk2, 0, kk2+kk2+nn2-kk2, nn2);
349 appendVector(hh2, hh2, h7, nn2, nn2);
351 appendVector(DD2, DD2, D7, kk2, kk2);
352 appendBlockMatrix(JJ2, JJ2, J7, kk2, kk2, kk2, kk2);
355 int n8 = nn2+nn3; int k8 = kk2+kk3; int **G8; int **GBar8; int *h8; int Q8; int *D8; int **J8;
356 G8 = calloc(nn2+nn3,sizeof(int*));
357 GBar8 = calloc(nn2+nn3,sizeof(int*));
358 h8 = calloc(nn2+nn3,sizeof(int));
359 D8 = calloc(kk2+kk3,sizeof(int));
360 J8 = calloc(kk2+kk3,sizeof(int*));
361 for(i=0; i<nn2+nn3; i++) {
362 G8[i] = calloc(nn2+nn3, sizeof(int)); GBar8[i] = calloc(nn2+nn3, sizeof(int));
364 for(i=0; i<kk2+kk3; i++)
365 J8[i] = calloc(kk2+kk3, sizeof(int));
366 //appendBlockMatrix(GG3, GG2, G10, nn3, nn3, nn2, nn2);
367 addSubMatrix(GG3, G8, kk3, nn3, 0, 0, 0, 0);
368 addSubMatrix(GG2, G8, kk2, nn2, 0, 0, kk3, nn3);
369 addSubMatrix(GG3, G8, nn3-kk3, nn3, kk3, 0, kk3+kk2, 0); // put the out-of-subspace vectors at the end
370 addSubMatrix(GG2, G8, nn2-kk2, nn2, kk2, 0, kk3+kk2+nn3-kk3, nn2);
371 //appendBlockMatrix(GGBar3, GGBar2, GBar10, nn3, nn3, nn2, nn2);
372 addSubMatrix(GGBar3, GBar8, kk3, nn3, 0, 0, 0, 0);
373 addSubMatrix(GGBar2, GBar8, kk2, nn2, 0, 0, kk3, nn3);
374 addSubMatrix(GGBar3, GBar8, nn3-kk3, nn3, kk3, 0, kk3+kk2, 0); // put the out-of-subspace vectors at the end
375 addSubMatrix(GGBar2, GBar8, nn2-kk2, nn2, kk2, 0, kk3+kk2+nn3-kk3, nn2);
376 appendVector(hh3, hh2, h8, nn3, nn2);
378 appendVector(DD3, DD2, D8, kk3, kk2);
379 appendBlockMatrix(JJ3, JJ2, J8, kk3, kk3, kk2, kk2);
382 int n9 = nn2+nn4; int k9 = kk2+kk4; int **G9; int **GBar9; int *h9; int Q9; int *D9; int **J9;
383 G9 = calloc(nn2+nn4,sizeof(int*));
384 GBar9 = calloc(nn2+nn4,sizeof(int*));
385 h9 = calloc(nn2+nn4,sizeof(int));
386 D9 = calloc(kk2+kk4,sizeof(int));
387 J9 = calloc(kk2+kk4,sizeof(int*));
388 for(i=0; i<nn2+nn4; i++) {
389 G9[i] = calloc(nn2+nn4, sizeof(int)); GBar9[i] = calloc(nn2+nn4, sizeof(int));
391 for(i=0; i<kk2+kk4; i++)
392 J9[i] = calloc(kk2+kk4, sizeof(int));
393 //appendBlockMatrix(GG4, GG2, G11, nn4, nn4, nn2, nn2);
394 addSubMatrix(GG4, G9, kk4, nn4, 0, 0, 0, 0);
395 addSubMatrix(GG2, G9, kk2, nn2, 0, 0, kk4, nn4);
396 addSubMatrix(GG4, G9, nn4-kk4, nn4, kk4, 0, kk4+kk2, 0); // put the out-of-subspace vectors at the end
397 addSubMatrix(GG2, G9, nn2-kk2, nn2, kk2, 0, kk4+kk2+nn4-kk4, nn2);
398 //appendBlockMatrix(GGBar4, GGBar2, GBar11, nn4, nn4, nn2, nn2);
399 addSubMatrix(GGBar4, GBar9, kk4, nn4, 0, 0, 0, 0);
400 addSubMatrix(GGBar2, GBar9, kk2, nn2, 0, 0, kk4, nn4);
401 addSubMatrix(GGBar4, GBar9, nn4-kk4, nn4, kk4, 0, kk4+kk2, 0); // put the out-of-subspace vectors at the end
402 addSubMatrix(GGBar2, GBar9, nn2-kk2, nn2, kk2, 0, kk4+kk2+nn4-kk4, nn2);
403 appendVector(hh4, hh2, h9, nn4, nn2);
405 appendVector(DD4, DD2, D9, kk4, kk2);
406 appendBlockMatrix(JJ4, JJ2, J9, kk4, kk4, kk2, kk2);
409 int n10 = nn2+nn5; int k10 = kk2+kk5; int **G10; int **GBar10; int *h10; int Q10; int *D10; int **J10;
410 G10 = calloc(nn2+nn5,sizeof(int*));
411 GBar10 = calloc(nn2+nn5,sizeof(int*));
412 h10 = calloc(nn2+nn5,sizeof(int));
413 D10 = calloc(kk2+kk5,sizeof(int));
414 J10 = calloc(kk2+kk5,sizeof(int*));
415 for(i=0; i<nn2+nn5; i++) {
416 G10[i] = calloc(nn2+nn5, sizeof(int)); GBar10[i] = calloc(nn2+nn5, sizeof(int));
418 for(i=0; i<kk2+kk5; i++)
419 J10[i] = calloc(kk2+kk5, sizeof(int));
420 //appendBlockMatrix(GG5, GG2, G12, nn5, nn5, nn2, nn2);
421 addSubMatrix(GG5, G10, kk5, nn5, 0, 0, 0, 0);
422 addSubMatrix(GG2, G10, kk2, nn2, 0, 0, kk5, nn5);
423 addSubMatrix(GG5, G10, nn5-kk5, nn5, kk5, 0, kk5+kk2, 0); // put the out-of-subspace vectors at the end
424 addSubMatrix(GG2, G10, nn2-kk2, nn2, kk2, 0, kk5+kk2+nn5-kk5, nn2);
425 //appendBlockMatrix(GGBar5, GGBar2, GBar12, nn5, nn5, nn2, nn2);
426 addSubMatrix(GGBar5, GBar10, kk5, nn5, 0, 0, 0, 0);
427 addSubMatrix(GGBar2, GBar10, kk2, nn2, 0, 0, kk5, nn5);
428 addSubMatrix(GGBar5, GBar10, nn5-kk5, nn5, kk5, 0, kk5+kk2, 0); // put the out-of-subspace vectors at the end
429 addSubMatrix(GGBar2, GBar10, nn2-kk2, nn2, kk2, 0, kk5+kk2+nn5-kk5, nn2);
430 appendVector(hh5, hh2, h10, nn5, nn2);
432 appendVector(DD5, DD2, D10, kk5, kk2);
433 appendBlockMatrix(JJ5, JJ2, J10, kk5, kk5, kk2, kk2);
436 int n11 = nn2+nn6; int k11 = kk2+kk6; int **G11; int **GBar11; int *h11; int Q11; int *D11; int **J11;
437 G11 = calloc(nn2+nn6,sizeof(int*));
438 GBar11 = calloc(nn2+nn6,sizeof(int*));
439 h11 = calloc(nn2+nn6,sizeof(int));
440 D11 = calloc(kk2+kk6,sizeof(int));
441 J11 = calloc(kk2+kk6,sizeof(int*));
442 for(i=0; i<nn2+nn6; i++) {
443 G11[i] = calloc(nn2+nn6, sizeof(int)); GBar11[i] = calloc(nn2+nn6, sizeof(int));
445 for(i=0; i<kk2+kk6; i++)
446 J11[i] = calloc(kk2+kk6, sizeof(int));
447 //appendBlockMatrix(GG6, GG2, G13, nn6, nn6, nn2, nn2);
448 addSubMatrix(GG6, G11, kk6, nn6, 0, 0, 0, 0);
449 addSubMatrix(GG2, G11, kk2, nn2, 0, 0, kk6, nn6);
450 addSubMatrix(GG6, G11, nn6-kk6, nn6, kk6, 0, kk6+kk2, 0); // put the out-of-subspace vectors at the end
451 addSubMatrix(GG2, G11, nn2-kk2, nn2, kk2, 0, kk6+kk2+nn6-kk6, nn2);
452 //appendBlockMatrix(GGBar6, GGBar2, GBar13, nn6, nn6, nn2, nn2);
453 addSubMatrix(GGBar6, GBar11, kk6, nn6, 0, 0, 0, 0);
454 addSubMatrix(GGBar2, GBar11, kk2, nn2, 0, 0, kk6, nn6);
455 addSubMatrix(GGBar6, GBar11, nn6-kk6, nn6, kk6, 0, kk6+kk2, 0); // put the out-of-subspace vectors at the end
456 addSubMatrix(GGBar2, GBar11, nn2-kk2, nn2, kk2, 0, kk6+kk2+nn6-kk6, nn2);
457 appendVector(hh6, hh2, h11, nn6, nn2);
459 appendVector(DD6, DD2, D11, kk6, kk2);
460 appendBlockMatrix(JJ6, JJ2, J11, kk6, kk6, kk2, kk2);
463 int n12 = nn2+nn7; int k12 = kk2+kk7; int **G12; int **GBar12; int *h12; int Q12; int *D12; int **J12;
464 G12 = calloc(nn2+nn7,sizeof(int*));
465 GBar12 = calloc(nn2+nn7,sizeof(int*));
466 h12 = calloc(nn2+nn7,sizeof(int));
467 D12 = calloc(kk2+kk7,sizeof(int));
468 J12 = calloc(kk2+kk7,sizeof(int*));
469 for(i=0; i<nn2+nn7; i++) {
470 G12[i] = calloc(nn2+nn7, sizeof(int)); GBar12[i] = calloc(nn2+nn7, sizeof(int));
472 for(i=0; i<kk2+kk7; i++)
473 J12[i] = calloc(kk2+kk7, sizeof(int));
474 //appendBlockMatrix(GG7, GG2, G14, nn7, nn7, nn2, nn2);
475 addSubMatrix(GG7, G12, kk7, nn7, 0, 0, 0, 0);
476 addSubMatrix(GG2, G12, kk2, nn2, 0, 0, kk7, nn7);
477 addSubMatrix(GG7, G12, nn7-kk7, nn7, kk7, 0, kk7+kk2, 0); // put the out-of-subspace vectors at the end
478 addSubMatrix(GG2, G12, nn2-kk2, nn2, kk2, 0, kk7+kk2+nn7-kk7, nn2);
479 //appendBlockMatrix(GGBar7, GGBar2, GBar14, nn7, nn7, nn2, nn2);
480 addSubMatrix(GGBar7, GBar12, kk7, nn7, 0, 0, 0, 0);
481 addSubMatrix(GGBar2, GBar12, kk2, nn2, 0, 0, kk7, nn7);
482 addSubMatrix(GGBar7, GBar12, nn7-kk7, nn7, kk7, 0, kk7+kk2, 0); // put the out-of-subspace vectors at the end
483 addSubMatrix(GGBar2, GBar12, nn2-kk2, nn2, kk2, 0, kk7+kk2+nn7-kk7, nn2);
484 appendVector(hh7, hh2, h12, nn7, nn2);
486 appendVector(DD7, DD2, D12, kk7, kk2);
487 appendBlockMatrix(JJ7, JJ2, J12, kk7, kk7, kk2, kk2);
490 int n13 = nn3+nn1; int k13 = kk3+kk1; int **G13; int **GBar13; int *h13; int Q13; int *D13; int **J13;
491 G13 = calloc(nn3+nn1,sizeof(int*));
492 GBar13 = calloc(nn3+nn1,sizeof(int*));
493 h13 = calloc(nn3+nn1,sizeof(int));
494 D13 = calloc(kk3+kk1,sizeof(int));
495 J13 = calloc(kk3+kk1,sizeof(int*));
496 for(i=0; i<nn3+nn1; i++) {
497 G13[i] = calloc(nn3+nn1, sizeof(int)); GBar13[i] = calloc(nn3+nn1, sizeof(int));
499 for(i=0; i<kk3+kk1; i++)
500 J13[i] = calloc(kk3+kk1, sizeof(int));
501 //appendBlockMatrix(GG1, GG3, G15, nn1, nn1, nn3, nn3);
502 addSubMatrix(GG1, G13, kk1, nn1, 0, 0, 0, 0);
503 addSubMatrix(GG3, G13, kk3, nn3, 0, 0, kk1, nn1);
504 addSubMatrix(GG1, G13, nn1-kk1, nn1, kk1, 0, kk1+kk3, 0); // put the out-of-subspace vectors at the end
505 addSubMatrix(GG3, G13, nn3-kk3, nn3, kk3, 0, kk1+kk3+nn1-kk1, nn3);
506 //appendBlockMatrix(GGBar1, GGBar3, GBar15, nn1, nn1, nn3, nn3);
507 addSubMatrix(GGBar1, GBar13, kk1, nn1, 0, 0, 0, 0);
508 addSubMatrix(GGBar3, GBar13, kk3, nn3, 0, 0, kk1, nn1);
509 addSubMatrix(GGBar1, GBar13, nn1-kk1, nn1, kk1, 0, kk1+kk3, 0); // put the out-of-subspace vectors at the end
510 addSubMatrix(GGBar3, GBar13, nn3-kk3, nn3, kk3, 0, kk1+kk3+nn1-kk1, nn3);
511 appendVector(hh1, hh3, h13, nn1, nn3);
513 appendVector(DD1, DD3, D13, kk1, kk3);
514 appendBlockMatrix(JJ1, JJ3, J13, kk1, kk1, kk3, kk3);
517 int n14 = nn3+nn2; int k14 = kk3+kk2; int **G14; int **GBar14; int *h14; int Q14; int *D14; int **J14;
518 G14 = calloc(nn3+nn2,sizeof(int*));
519 GBar14 = calloc(nn3+nn2,sizeof(int*));
520 h14 = calloc(nn3+nn2,sizeof(int));
521 D14 = calloc(kk3+kk2,sizeof(int));
522 J14 = calloc(kk3+kk2,sizeof(int*));
523 for(i=0; i<nn3+nn2; i++) {
524 G14[i] = calloc(nn3+nn2, sizeof(int)); GBar14[i] = calloc(nn3+nn2, sizeof(int));
526 for(i=0; i<kk3+kk2; i++)
527 J14[i] = calloc(kk3+kk2, sizeof(int));
528 //appendBlockMatrix(GG2, GG3, G16, nn2, nn2, nn3, nn3);
529 addSubMatrix(GG2, G14, kk2, nn2, 0, 0, 0, 0);
530 addSubMatrix(GG3, G14, kk3, nn3, 0, 0, kk2, nn2);
531 addSubMatrix(GG2, G14, nn2-kk2, nn2, kk2, 0, kk2+kk3, 0); // put the out-of-subspace vectors at the end
532 addSubMatrix(GG3, G14, nn3-kk3, nn3, kk3, 0, kk2+kk3+nn2-kk2, nn3);
533 //appendBlockMatrix(GGBar2, GGBar3, GBar16, nn2, nn2, nn3, nn3);
534 addSubMatrix(GGBar2, GBar14, kk2, nn2, 0, 0, 0, 0);
535 addSubMatrix(GGBar3, GBar14, kk3, nn3, 0, 0, kk2, nn2);
536 addSubMatrix(GGBar2, GBar14, nn2-kk2, nn2, kk2, 0, kk2+kk3, 0); // put the out-of-subspace vectors at the end
537 addSubMatrix(GGBar3, GBar14, nn3-kk3, nn3, kk3, 0, kk2+kk3+nn2-kk2, nn3);
538 appendVector(hh2, hh3, h14, nn2, nn3);
540 appendVector(DD2, DD3, D14, kk2, kk3);
541 appendBlockMatrix(JJ2, JJ3, J14, kk2, kk2, kk3, kk3);
544 int n15 = nn3+nn3; int k15 = kk3+kk3; int **G15; int **GBar15; int *h15; int Q15; int *D15; int **J15;
545 G15 = calloc(nn3+nn3,sizeof(int*));
546 GBar15 = calloc(nn3+nn3,sizeof(int*));
547 h15 = calloc(nn3+nn3,sizeof(int));
548 D15 = calloc(kk3+kk3,sizeof(int));
549 J15 = calloc(kk3+kk3,sizeof(int*));
550 for(i=0; i<nn3+nn3; i++) {
551 G15[i] = calloc(nn3+nn3, sizeof(int)); GBar15[i] = calloc(nn3+nn3, sizeof(int));
553 for(i=0; i<kk3+kk3; i++)
554 J15[i] = calloc(kk3+kk3, sizeof(int));
555 //appendBlockMatrix(GG3, GG3, G17, nn3, nn3, nn3, nn3);
556 addSubMatrix(GG3, G15, kk3, nn3, 0, 0, 0, 0);
557 addSubMatrix(GG3, G15, kk3, nn3, 0, 0, kk3, nn3);
558 addSubMatrix(GG3, G15, nn3-kk3, nn3, kk3, 0, kk3+kk3, 0); // put the out-of-subspace vectors at the end
559 addSubMatrix(GG3, G15, nn3-kk3, nn3, kk3, 0, kk3+kk3+nn3-kk3, nn3);
560 //appendBlockMatrix(GGBar3, GGBar3, GBar17, nn3, nn3, nn3, nn3);
561 addSubMatrix(GGBar3, GBar15, kk3, nn3, 0, 0, 0, 0);
562 addSubMatrix(GGBar3, GBar15, kk3, nn3, 0, 0, kk3, nn3);
563 addSubMatrix(GGBar3, GBar15, nn3-kk3, nn3, kk3, 0, kk3+kk3, 0); // put the out-of-subspace vectors at the end
564 addSubMatrix(GGBar3, GBar15, nn3-kk3, nn3, kk3, 0, kk3+kk3+nn3-kk3, nn3);
565 appendVector(hh3, hh3, h15, nn3, nn3);
567 appendVector(DD3, DD3, D15, kk3, kk3);
568 appendBlockMatrix(JJ3, JJ3, J15, kk3, kk3, kk3, kk3);
571 /*int n18 = nn3+nn4; int k18 = kk3+kk4; int **G18; int **GBar18; int *h18; int Q18; int *D18; int **J18;
572 G18 = calloc(nn3+nn4,sizeof(int*));
573 GBar18 = calloc(nn3+nn4,sizeof(int*));
574 h18 = calloc(nn3+nn4,sizeof(int));
575 D18 = calloc(kk3+kk4,sizeof(int));
576 J18 = calloc(kk3+kk4,sizeof(int*));
577 for(i=0; i<nn3+nn4; i++) {
578 G18[i] = calloc(nn3+nn4, sizeof(int)); GBar18[i] = calloc(nn3+nn4, sizeof(int));
580 for(i=0; i<kk3+kk4; i++)
581 J18[i] = calloc(kk3+kk4, sizeof(int));
582 //appendBlockMatrix(GG4, GG3, G18, nn4, nn4, nn3, nn3);
583 addSubMatrix(GG4, G18, kk4, nn4, 0, 0, 0, 0);
584 addSubMatrix(GG3, G18, kk3, nn3, 0, 0, kk4, nn4);
585 addSubMatrix(GG4, G18, nn4-kk4, nn4, kk4, 0, kk4+kk3, 0); // put the out-of-subspace vectors at the end
586 addSubMatrix(GG3, G18, nn3-kk3, nn3, kk3, 0, kk4+kk3+nn4-kk4, nn3);
587 //appendBlockMatrix(GGBar4, GGBar3, GBar18, nn4, nn4, nn3, nn3);
588 addSubMatrix(GGBar4, GBar18, kk4, nn4, 0, 0, 0, 0);
589 addSubMatrix(GGBar3, GBar18, kk3, nn3, 0, 0, kk4, nn4);
590 addSubMatrix(GGBar4, GBar18, nn4-kk4, nn4, kk4, 0, kk4+kk3, 0); // put the out-of-subspace vectors at the end
591 addSubMatrix(GGBar3, GBar18, nn3-kk3, nn3, kk3, 0, kk4+kk3+nn4-kk4, nn3);
592 appendVector(hh4, hh3, h18, nn4, nn3);
594 appendVector(DD4, DD3, D18, kk4, kk3);
595 appendBlockMatrix(JJ4, JJ3, J18, kk4, kk4, kk3, kk3);*/
598 int n16 = nn3+nn5; int k16 = kk3+kk5; int **G16; int **GBar16; int *h16; int Q16; int *D16; int **J16;
599 G16 = calloc(nn3+nn5,sizeof(int*));
600 GBar16 = calloc(nn3+nn5,sizeof(int*));
601 h16 = calloc(nn3+nn5,sizeof(int));
602 D16 = calloc(kk3+kk5,sizeof(int));
603 J16 = calloc(kk3+kk5,sizeof(int*));
604 for(i=0; i<nn3+nn5; i++) {
605 G16[i] = calloc(nn3+nn5, sizeof(int)); GBar16[i] = calloc(nn3+nn5, sizeof(int));
607 for(i=0; i<kk3+kk5; i++)
608 J16[i] = calloc(kk3+kk5, sizeof(int));
609 //appendBlockMatrix(GG5, GG3, G19, nn5, nn5, nn3, nn3);
610 addSubMatrix(GG5, G16, kk5, nn5, 0, 0, 0, 0);
611 addSubMatrix(GG3, G16, kk3, nn3, 0, 0, kk5, nn5);
612 addSubMatrix(GG5, G16, nn5-kk5, nn5, kk5, 0, kk5+kk3, 0); // put the out-of-subspace vectors at the end
613 addSubMatrix(GG3, G16, nn3-kk3, nn3, kk3, 0, kk5+kk3+nn5-kk5, nn3);
614 //appendBlockMatrix(GGBar5, GGBar3, GBar19, nn5, nn5, nn3, nn3);
615 addSubMatrix(GGBar5, GBar16, kk5, nn5, 0, 0, 0, 0);
616 addSubMatrix(GGBar3, GBar16, kk3, nn3, 0, 0, kk5, nn5);
617 addSubMatrix(GGBar5, GBar16, nn5-kk5, nn5, kk5, 0, kk5+kk3, 0); // put the out-of-subspace vectors at the end
618 addSubMatrix(GGBar3, GBar16, nn3-kk3, nn3, kk3, 0, kk5+kk3+nn5-kk5, nn3);
619 appendVector(hh5, hh3, h16, nn5, nn3);
621 appendVector(DD5, DD3, D16, kk5, kk3);
622 appendBlockMatrix(JJ5, JJ3, J16, kk5, kk5, kk3, kk3);
625 int n17 = nn3+nn6; int k17 = kk3+kk6; int **G17; int **GBar17; int *h17; int Q17; int *D17; int **J17;
626 G17 = calloc(nn3+nn6,sizeof(int*));
627 GBar17 = calloc(nn3+nn6,sizeof(int*));
628 h17 = calloc(nn3+nn6,sizeof(int));
629 D17 = calloc(kk3+kk6,sizeof(int));
630 J17 = calloc(kk3+kk6,sizeof(int*));
631 for(i=0; i<nn3+nn6; i++) {
632 G17[i] = calloc(nn3+nn6, sizeof(int)); GBar17[i] = calloc(nn3+nn6, sizeof(int));
634 for(i=0; i<kk3+kk6; i++)
635 J17[i] = calloc(kk3+kk6, sizeof(int));
636 //appendBlockMatrix(GG6, GG3, G20, nn6, nn6, nn3, nn3);
637 addSubMatrix(GG6, G17, kk6, nn6, 0, 0, 0, 0);
638 addSubMatrix(GG3, G17, kk3, nn3, 0, 0, kk6, nn6);
639 addSubMatrix(GG6, G17, nn6-kk6, nn6, kk6, 0, kk6+kk3, 0); // put the out-of-subspace vectors at the end
640 addSubMatrix(GG3, G17, nn3-kk3, nn3, kk3, 0, kk6+kk3+nn6-kk6, nn3);
641 //appendBlockMatrix(GGBar6, GGBar3, GBar20, nn6, nn6, nn3, nn3);
642 addSubMatrix(GGBar6, GBar17, kk6, nn6, 0, 0, 0, 0);
643 addSubMatrix(GGBar3, GBar17, kk3, nn3, 0, 0, kk6, nn6);
644 addSubMatrix(GGBar6, GBar17, nn6-kk6, nn6, kk6, 0, kk6+kk3, 0); // put the out-of-subspace vectors at the end
645 addSubMatrix(GGBar3, GBar17, nn3-kk3, nn3, kk3, 0, kk6+kk3+nn6-kk6, nn3);
646 appendVector(hh6, hh3, h17, nn6, nn3);
648 appendVector(DD6, DD3, D17, kk6, kk3);
649 appendBlockMatrix(JJ6, JJ3, J17, kk6, kk6, kk3, kk3);
652 int n18 = nn3+nn7; int k18 = kk3+kk7; int **G18; int **GBar18; int *h18; int Q18; int *D18; int **J18;
653 G18 = calloc(nn3+nn7,sizeof(int*));
654 GBar18 = calloc(nn3+nn7,sizeof(int*));
655 h18 = calloc(nn3+nn7,sizeof(int));
656 D18 = calloc(kk3+kk7,sizeof(int));
657 J18 = calloc(kk3+kk7,sizeof(int*));
658 for(i=0; i<nn3+nn7; i++) {
659 G18[i] = calloc(nn3+nn7, sizeof(int)); GBar18[i] = calloc(nn3+nn7, sizeof(int));
661 for(i=0; i<kk3+kk7; i++)
662 J18[i] = calloc(kk3+kk7, sizeof(int));
663 //appendBlockMatrix(GG7, GG3, G21, nn7, nn7, nn3, nn3);
664 addSubMatrix(GG7, G18, kk7, nn7, 0, 0, 0, 0);
665 addSubMatrix(GG3, G18, kk3, nn3, 0, 0, kk7, nn7);
666 addSubMatrix(GG7, G18, nn7-kk7, nn7, kk7, 0, kk7+kk3, 0); // put the out-of-subspace vectors at the end
667 addSubMatrix(GG3, G18, nn3-kk3, nn3, kk3, 0, kk7+kk3+nn7-kk7, nn3);
668 //appendBlockMatrix(GGBar7, GGBar3, GBar21, nn7, nn7, nn3, nn3);
669 addSubMatrix(GGBar7, GBar18, kk7, nn7, 0, 0, 0, 0);
670 addSubMatrix(GGBar3, GBar18, kk3, nn3, 0, 0, kk7, nn7);
671 addSubMatrix(GGBar7, GBar18, nn7-kk7, nn7, kk7, 0, kk7+kk3, 0); // put the out-of-subspace vectors at the end
672 addSubMatrix(GGBar3, GBar18, nn3-kk3, nn3, kk3, 0, kk7+kk3+nn7-kk7, nn3);
673 appendVector(hh7, hh3, h18, nn7, nn3);
675 appendVector(DD7, DD3, D18, kk7, kk3);
676 appendBlockMatrix(JJ7, JJ3, J18, kk7, kk7, kk3, kk3);
679 int n19 = nn4+nn1; int k19 = kk4+kk1; int **G19; int **GBar19; int *h19; int Q19; int *D19; int **J19;
680 G19 = calloc(nn4+nn1,sizeof(int*));
681 GBar19 = calloc(nn4+nn1,sizeof(int*));
682 h19 = calloc(nn4+nn1,sizeof(int));
683 D19 = calloc(kk4+kk1,sizeof(int));
684 J19 = calloc(kk4+kk1,sizeof(int*));
685 for(i=0; i<nn4+nn1; i++) {
686 G19[i] = calloc(nn4+nn1, sizeof(int)); GBar19[i] = calloc(nn4+nn1, sizeof(int));
688 for(i=0; i<kk4+kk1; i++)
689 J19[i] = calloc(kk4+kk1, sizeof(int));
690 //appendBlockMatrix(GG1, GG4, G22, nn1, nn1, nn4, nn4);
691 addSubMatrix(GG1, G19, kk1, nn1, 0, 0, 0, 0);
692 addSubMatrix(GG4, G19, kk4, nn4, 0, 0, kk1, nn1);
693 addSubMatrix(GG1, G19, nn1-kk1, nn1, kk1, 0, kk1+kk4, 0); // put the out-of-subspace vectors at the end
694 addSubMatrix(GG4, G19, nn4-kk4, nn4, kk4, 0, kk1+kk4+nn1-kk1, nn4);
695 //appendBlockMatrix(GGBar1, GGBar4, GBar22, nn1, nn1, nn4, nn4);
696 addSubMatrix(GGBar1, GBar19, kk1, nn1, 0, 0, 0, 0);
697 addSubMatrix(GGBar4, GBar19, kk4, nn4, 0, 0, kk1, nn1);
698 addSubMatrix(GGBar1, GBar19, nn1-kk1, nn1, kk1, 0, kk1+kk4, 0); // put the out-of-subspace vectors at the end
699 addSubMatrix(GGBar4, GBar19, nn4-kk4, nn4, kk4, 0, kk1+kk4+nn1-kk1, nn4);
700 appendVector(hh1, hh4, h19, nn1, nn4);
702 appendVector(DD1, DD4, D19, kk1, kk4);
703 appendBlockMatrix(JJ1, JJ4, J19, kk1, kk1, kk4, kk4);
706 int n20 = nn4+nn2; int k20 = kk4+kk2; int **G20; int **GBar20; int *h20; int Q20; int *D20; int **J20;
707 G20 = calloc(nn4+nn2,sizeof(int*));
708 GBar20 = calloc(nn4+nn2,sizeof(int*));
709 h20 = calloc(nn4+nn2,sizeof(int));
710 D20 = calloc(kk4+kk2,sizeof(int));
711 J20 = calloc(kk4+kk2,sizeof(int*));
712 for(i=0; i<nn4+nn2; i++) {
713 G20[i] = calloc(nn4+nn2, sizeof(int)); GBar20[i] = calloc(nn4+nn2, sizeof(int));
715 for(i=0; i<kk4+kk2; i++)
716 J20[i] = calloc(kk4+kk2, sizeof(int));
717 //appendBlockMatrix(GG2, GG4, G23, nn2, nn2, nn4, nn4);
718 addSubMatrix(GG2, G20, kk2, nn2, 0, 0, 0, 0);
719 addSubMatrix(GG4, G20, kk4, nn4, 0, 0, kk2, nn2);
720 addSubMatrix(GG2, G20, nn2-kk2, nn2, kk2, 0, kk2+kk4, 0); // put the out-of-subspace vectors at the end
721 addSubMatrix(GG4, G20, nn4-kk4, nn4, kk4, 0, kk2+kk4+nn2-kk2, nn4);
722 //appendBlockMatrix(GGBar2, GGBar4, GBar23, nn2, nn2, nn4, nn4);
723 addSubMatrix(GGBar2, GBar20, kk2, nn2, 0, 0, 0, 0);
724 addSubMatrix(GGBar4, GBar20, kk4, nn4, 0, 0, kk2, nn2);
725 addSubMatrix(GGBar2, GBar20, nn2-kk2, nn2, kk2, 0, kk2+kk4, 0); // put the out-of-subspace vectors at the end
726 addSubMatrix(GGBar4, GBar20, nn4-kk4, nn4, kk4, 0, kk2+kk4+nn2-kk2, nn4);
727 appendVector(hh2, hh4, h20, nn2, nn4);
729 appendVector(DD2, DD4, D20, kk2, kk4);
730 appendBlockMatrix(JJ2, JJ4, J20, kk2, kk2, kk4, kk4);
733 /*int n24 = nn4+nn3; int k24 = kk4+kk3; int **G24; int **GBar24; int *h24; int Q24; int *D24; int **J24;
734 G24 = calloc(nn4+nn3,sizeof(int*));
735 GBar24 = calloc(nn4+nn3,sizeof(int*));
736 h24 = calloc(nn4+nn3,sizeof(int));
737 D24 = calloc(kk4+kk3,sizeof(int));
738 J24 = calloc(kk4+kk3,sizeof(int*));
739 for(i=0; i<nn4+nn3; i++) {
740 G24[i] = calloc(nn4+nn3, sizeof(int)); GBar24[i] = calloc(nn4+nn3, sizeof(int));
742 for(i=0; i<kk4+kk3; i++)
743 J24[i] = calloc(kk4+kk3, sizeof(int));
744 //appendBlockMatrix(GG3, GG4, G24, nn3, nn3, nn4, nn4);
745 addSubMatrix(GG3, G24, kk3, nn3, 0, 0, 0, 0);
746 addSubMatrix(GG4, G24, kk4, nn4, 0, 0, kk3, nn3);
747 addSubMatrix(GG3, G24, nn3-kk3, nn3, kk3, 0, kk3+kk4, 0); // put the out-of-subspace vectors at the end
748 addSubMatrix(GG4, G24, nn4-kk4, nn4, kk4, 0, kk3+kk4+nn3-kk3, nn4);
749 //appendBlockMatrix(GGBar3, GGBar4, GBar24, nn3, nn3, nn4, nn4);
750 addSubMatrix(GGBar3, GBar24, kk3, nn3, 0, 0, 0, 0);
751 addSubMatrix(GGBar4, GBar24, kk4, nn4, 0, 0, kk3, nn3);
752 addSubMatrix(GGBar3, GBar24, nn3-kk3, nn3, kk3, 0, kk3+kk4, 0); // put the out-of-subspace vectors at the end
753 addSubMatrix(GGBar4, GBar24, nn4-kk4, nn4, kk4, 0, kk3+kk4+nn3-kk3, nn4);
754 appendVector(hh3, hh4, h24, nn3, nn4);
756 appendVector(DD3, DD4, D24, kk3, kk4);
757 appendBlockMatrix(JJ3, JJ4, J24, kk3, kk3, kk4, kk4);*/
760 int n21 = nn4+nn4; int k21 = kk4+kk4; int **G21; int **GBar21; int *h21; int Q21; int *D21; int **J21;
761 G21 = calloc(nn4+nn4,sizeof(int*));
762 GBar21 = calloc(nn4+nn4,sizeof(int*));
763 h21 = calloc(nn4+nn4,sizeof(int));
764 D21 = calloc(kk4+kk4,sizeof(int));
765 J21 = calloc(kk4+kk4,sizeof(int*));
766 for(i=0; i<nn4+nn4; i++) {
767 G21[i] = calloc(nn4+nn4, sizeof(int)); GBar21[i] = calloc(nn4+nn4, sizeof(int));
769 for(i=0; i<kk4+kk4; i++)
770 J21[i] = calloc(kk4+kk4, sizeof(int));
771 //appendBlockMatrix(GG4, GG4, G25, nn4, nn4, nn4, nn4);
772 addSubMatrix(GG4, G21, kk4, nn4, 0, 0, 0, 0);
773 addSubMatrix(GG4, G21, kk4, nn4, 0, 0, kk4, nn4);
774 addSubMatrix(GG4, G21, nn4-kk4, nn4, kk4, 0, kk4+kk4, 0); // put the out-of-subspace vectors at the end
775 addSubMatrix(GG4, G21, nn4-kk4, nn4, kk4, 0, kk4+kk4+nn4-kk4, nn4);
776 //appendBlockMatrix(GGBar4, GGBar4, GBar25, nn4, nn4, nn4, nn4);
777 addSubMatrix(GGBar4, GBar21, kk4, nn4, 0, 0, 0, 0);
778 addSubMatrix(GGBar4, GBar21, kk4, nn4, 0, 0, kk4, nn4);
779 addSubMatrix(GGBar4, GBar21, nn4-kk4, nn4, kk4, 0, kk4+kk4, 0); // put the out-of-subspace vectors at the end
780 addSubMatrix(GGBar4, GBar21, nn4-kk4, nn4, kk4, 0, kk4+kk4+nn4-kk4, nn4);
781 appendVector(hh4, hh4, h21, nn4, nn4);
783 appendVector(DD4, DD4, D21, kk4, kk4);
784 appendBlockMatrix(JJ4, JJ4, J21, kk4, kk4, kk4, kk4);
787 int n22 = nn4+nn5; int k22 = kk4+kk5; int **G22; int **GBar22; int *h22; int Q22; int *D22; int **J22;
788 G22 = calloc(nn4+nn5,sizeof(int*));
789 GBar22 = calloc(nn4+nn5,sizeof(int*));
790 h22 = calloc(nn4+nn5,sizeof(int));
791 D22 = calloc(kk4+kk5,sizeof(int));
792 J22 = calloc(kk4+kk5,sizeof(int*));
793 for(i=0; i<nn4+nn5; i++) {
794 G22[i] = calloc(nn4+nn5, sizeof(int)); GBar22[i] = calloc(nn4+nn5, sizeof(int));
796 for(i=0; i<kk4+kk5; i++)
797 J22[i] = calloc(kk4+kk5, sizeof(int));
798 //appendBlockMatrix(GG5, GG4, G26, nn5, nn5, nn4, nn4);
799 addSubMatrix(GG5, G22, kk5, nn5, 0, 0, 0, 0);
800 addSubMatrix(GG4, G22, kk4, nn4, 0, 0, kk5, nn5);
801 addSubMatrix(GG5, G22, nn5-kk5, nn5, kk5, 0, kk5+kk4, 0); // put the out-of-subspace vectors at the end
802 addSubMatrix(GG4, G22, nn4-kk4, nn4, kk4, 0, kk5+kk4+nn5-kk5, nn4);
803 //appendBlockMatrix(GGBar5, GGBar4, GBar26, nn5, nn5, nn4, nn4);
804 addSubMatrix(GGBar5, GBar22, kk5, nn5, 0, 0, 0, 0);
805 addSubMatrix(GGBar4, GBar22, kk4, nn4, 0, 0, kk5, nn5);
806 addSubMatrix(GGBar5, GBar22, nn5-kk5, nn5, kk5, 0, kk5+kk4, 0); // put the out-of-subspace vectors at the end
807 addSubMatrix(GGBar4, GBar22, nn4-kk4, nn4, kk4, 0, kk5+kk4+nn5-kk5, nn4);
808 appendVector(hh5, hh4, h22, nn5, nn4);
810 appendVector(DD5, DD4, D22, kk5, kk4);
811 appendBlockMatrix(JJ5, JJ4, J22, kk5, kk5, kk4, kk4);
814 int n23 = nn4+nn6; int k23 = kk4+kk6; int **G23; int **GBar23; int *h23; int Q23; int *D23; int **J23;
815 G23 = calloc(nn4+nn6,sizeof(int*));
816 GBar23 = calloc(nn4+nn6,sizeof(int*));
817 h23 = calloc(nn4+nn6,sizeof(int));
818 D23 = calloc(kk4+kk6,sizeof(int));
819 J23 = calloc(kk4+kk6,sizeof(int*));
820 for(i=0; i<nn4+nn6; i++) {
821 G23[i] = calloc(nn4+nn6, sizeof(int)); GBar23[i] = calloc(nn4+nn6, sizeof(int));
823 for(i=0; i<kk4+kk6; i++)
824 J23[i] = calloc(kk4+kk6, sizeof(int));
825 //appendBlockMatrix(GG6, GG4, G27, nn6, nn6, nn4, nn4);
826 addSubMatrix(GG6, G23, kk6, nn6, 0, 0, 0, 0);
827 addSubMatrix(GG4, G23, kk4, nn4, 0, 0, kk6, nn6);
828 addSubMatrix(GG6, G23, nn6-kk6, nn6, kk6, 0, kk6+kk4, 0); // put the out-of-subspace vectors at the end
829 addSubMatrix(GG4, G23, nn4-kk4, nn4, kk4, 0, kk6+kk4+nn6-kk6, nn4);
830 //appendBlockMatrix(GGBar6, GGBar4, GBar27, nn6, nn6, nn4, nn4);
831 addSubMatrix(GGBar6, GBar23, kk6, nn6, 0, 0, 0, 0);
832 addSubMatrix(GGBar4, GBar23, kk4, nn4, 0, 0, kk6, nn6);
833 addSubMatrix(GGBar6, GBar23, nn6-kk6, nn6, kk6, 0, kk6+kk4, 0); // put the out-of-subspace vectors at the end
834 addSubMatrix(GGBar4, GBar23, nn4-kk4, nn4, kk4, 0, kk6+kk4+nn6-kk6, nn4);
835 appendVector(hh6, hh4, h23, nn6, nn4);
837 appendVector(DD6, DD4, D23, kk6, kk4);
838 appendBlockMatrix(JJ6, JJ4, J23, kk6, kk6, kk4, kk4);
841 int n24 = nn4+nn7; int k24 = kk4+kk7; int **G24; int **GBar24; int *h24; int Q24; int *D24; int **J24;
842 G24 = calloc(nn4+nn7,sizeof(int*));
843 GBar24 = calloc(nn4+nn7,sizeof(int*));
844 h24 = calloc(nn4+nn7,sizeof(int));
845 D24 = calloc(kk4+kk7,sizeof(int));
846 J24 = calloc(kk4+kk7,sizeof(int*));
847 for(i=0; i<nn4+nn7; i++) {
848 G24[i] = calloc(nn4+nn7, sizeof(int)); GBar24[i] = calloc(nn4+nn7, sizeof(int));
850 for(i=0; i<kk4+kk7; i++)
851 J24[i] = calloc(kk4+kk7, sizeof(int));
852 //appendBlockMatrix(GG7, GG4, G28, nn7, nn7, nn4, nn4);
853 addSubMatrix(GG7, G24, kk7, nn7, 0, 0, 0, 0);
854 addSubMatrix(GG4, G24, kk4, nn4, 0, 0, kk7, nn7);
855 addSubMatrix(GG7, G24, nn7-kk7, nn7, kk7, 0, kk7+kk4, 0); // put the out-of-subspace vectors at the end
856 addSubMatrix(GG4, G24, nn4-kk4, nn4, kk4, 0, kk7+kk4+nn7-kk7, nn4);
857 //appendBlockMatrix(GGBar7, GGBar4, GBar28, nn7, nn7, nn4, nn4);
858 addSubMatrix(GGBar7, GBar24, kk7, nn7, 0, 0, 0, 0);
859 addSubMatrix(GGBar4, GBar24, kk4, nn4, 0, 0, kk7, nn7);
860 addSubMatrix(GGBar7, GBar24, nn7-kk7, nn7, kk7, 0, kk7+kk4, 0); // put the out-of-subspace vectors at the end
861 addSubMatrix(GGBar4, GBar24, nn4-kk4, nn4, kk4, 0, kk7+kk4+nn7-kk7, nn4);
862 appendVector(hh7, hh4, h24, nn7, nn4);
864 appendVector(DD7, DD4, D24, kk7, kk4);
865 appendBlockMatrix(JJ7, JJ4, J24, kk7, kk7, kk4, kk4);
868 int n25 = nn5+nn1; int k25 = kk5+kk1; int **G25; int **GBar25; int *h25; int Q25; int *D25; int **J25;
869 G25 = calloc(nn5+nn1,sizeof(int*));
870 GBar25 = calloc(nn5+nn1,sizeof(int*));
871 h25 = calloc(nn5+nn1,sizeof(int));
872 D25 = calloc(kk5+kk1,sizeof(int));
873 J25 = calloc(kk5+kk1,sizeof(int*));
874 for(i=0; i<nn5+nn1; i++) {
875 G25[i] = calloc(nn5+nn1, sizeof(int)); GBar25[i] = calloc(nn5+nn1, sizeof(int));
877 for(i=0; i<kk5+kk1; i++)
878 J25[i] = calloc(kk5+kk1, sizeof(int));
879 //appendBlockMatrix(GG1, GG5, G29, nn1, nn1, nn5, nn5);
880 addSubMatrix(GG1, G25, kk1, nn1, 0, 0, 0, 0);
881 addSubMatrix(GG5, G25, kk5, nn5, 0, 0, kk1, nn1);
882 addSubMatrix(GG1, G25, nn1-kk1, nn1, kk1, 0, kk1+kk5, 0); // put the out-of-subspace vectors at the end
883 addSubMatrix(GG5, G25, nn5-kk5, nn5, kk5, 0, kk1+kk5+nn1-kk1, nn5);
884 //appendBlockMatrix(GGBar1, GGBar5, GBar29, nn1, nn1, nn5, nn5);
885 addSubMatrix(GGBar1, GBar25, kk1, nn1, 0, 0, 0, 0);
886 addSubMatrix(GGBar5, GBar25, kk5, nn5, 0, 0, kk1, nn1);
887 addSubMatrix(GGBar1, GBar25, nn1-kk1, nn1, kk1, 0, kk1+kk5, 0); // put the out-of-subspace vectors at the end
888 addSubMatrix(GGBar5, GBar25, nn5-kk5, nn5, kk5, 0, kk1+kk5+nn1-kk1, nn5);
889 appendVector(hh1, hh5, h25, nn1, nn5);
891 appendVector(DD1, DD5, D25, kk1, kk5);
892 appendBlockMatrix(JJ1, JJ5, J25, kk1, kk1, kk5, kk5);
895 int n26 = nn5+nn2; int k26 = kk5+kk2; int **G26; int **GBar26; int *h26; int Q26; int *D26; int **J26;
896 G26 = calloc(nn5+nn2,sizeof(int*));
897 GBar26 = calloc(nn5+nn2,sizeof(int*));
898 h26 = calloc(nn5+nn2,sizeof(int));
899 D26 = calloc(kk5+kk2,sizeof(int));
900 J26 = calloc(kk5+kk2,sizeof(int*));
901 for(i=0; i<nn5+nn2; i++) {
902 G26[i] = calloc(nn5+nn2, sizeof(int)); GBar26[i] = calloc(nn5+nn2, sizeof(int));
904 for(i=0; i<kk5+kk2; i++)
905 J26[i] = calloc(kk5+kk2, sizeof(int));
906 //appendBlockMatrix(GG2, GG5, G30, nn2, nn2, nn5, nn5);
907 addSubMatrix(GG2, G26, kk2, nn2, 0, 0, 0, 0);
908 addSubMatrix(GG5, G26, kk5, nn5, 0, 0, kk2, nn2);
909 addSubMatrix(GG2, G26, nn2-kk2, nn2, kk2, 0, kk2+kk5, 0); // put the out-of-subspace vectors at the end
910 addSubMatrix(GG5, G26, nn5-kk5, nn5, kk5, 0, kk2+kk5+nn2-kk2, nn5);
911 //appendBlockMatrix(GGBar2, GGBar5, GBar30, nn2, nn2, nn5, nn5);
912 addSubMatrix(GGBar2, GBar26, kk2, nn2, 0, 0, 0, 0);
913 addSubMatrix(GGBar5, GBar26, kk5, nn5, 0, 0, kk2, nn2);
914 addSubMatrix(GGBar2, GBar26, nn2-kk2, nn2, kk2, 0, kk2+kk5, 0); // put the out-of-subspace vectors at the end
915 addSubMatrix(GGBar5, GBar26, nn5-kk5, nn5, kk5, 0, kk2+kk5+nn2-kk2, nn5);
916 appendVector(hh2, hh5, h26, nn2, nn5);
918 appendVector(DD2, DD5, D26, kk2, kk5);
919 appendBlockMatrix(JJ2, JJ5, J26, kk2, kk2, kk5, kk5);
922 int n27 = nn5+nn3; int k27 = kk5+kk3; int **G27; int **GBar27; int *h27; int Q27; int *D27; int **J27;
923 G27 = calloc(nn5+nn3,sizeof(int*));
924 GBar27 = calloc(nn5+nn3,sizeof(int*));
925 h27 = calloc(nn5+nn3,sizeof(int));
926 D27 = calloc(kk5+kk3,sizeof(int));
927 J27 = calloc(kk5+kk3,sizeof(int*));
928 for(i=0; i<nn5+nn3; i++) {
929 G27[i] = calloc(nn5+nn3, sizeof(int)); GBar27[i] = calloc(nn5+nn3, sizeof(int));
931 for(i=0; i<kk5+kk3; i++)
932 J27[i] = calloc(kk5+kk3, sizeof(int));
933 //appendBlockMatrix(GG3, GG5, G31, nn3, nn3, nn5, nn5);
934 addSubMatrix(GG3, G27, kk3, nn3, 0, 0, 0, 0);
935 addSubMatrix(GG5, G27, kk5, nn5, 0, 0, kk3, nn3);
936 addSubMatrix(GG3, G27, nn3-kk3, nn3, kk3, 0, kk3+kk5, 0); // put the out-of-subspace vectors at the end
937 addSubMatrix(GG5, G27, nn5-kk5, nn5, kk5, 0, kk3+kk5+nn3-kk3, nn5);
938 //appendBlockMatrix(GGBar3, GGBar5, GBar31, nn3, nn3, nn5, nn5);
939 addSubMatrix(GGBar3, GBar27, kk3, nn3, 0, 0, 0, 0);
940 addSubMatrix(GGBar5, GBar27, kk5, nn5, 0, 0, kk3, nn3);
941 addSubMatrix(GGBar3, GBar27, nn3-kk3, nn3, kk3, 0, kk3+kk5, 0); // put the out-of-subspace vectors at the end
942 addSubMatrix(GGBar5, GBar27, nn5-kk5, nn5, kk5, 0, kk3+kk5+nn3-kk3, nn5);
943 appendVector(hh3, hh5, h27, nn3, nn5);
945 appendVector(DD3, DD5, D27, kk3, kk5);
946 appendBlockMatrix(JJ3, JJ5, J27, kk3, kk3, kk5, kk5);
949 int n28 = nn5+nn4; int k28 = kk5+kk4; int **G28; int **GBar28; int *h28; int Q28; int *D28; int **J28;
950 G28 = calloc(nn5+nn4,sizeof(int*));
951 GBar28 = calloc(nn5+nn4,sizeof(int*));
952 h28 = calloc(nn5+nn4,sizeof(int));
953 D28 = calloc(kk5+kk4,sizeof(int));
954 J28 = calloc(kk5+kk4,sizeof(int*));
955 for(i=0; i<nn5+nn4; i++) {
956 G28[i] = calloc(nn5+nn4, sizeof(int)); GBar28[i] = calloc(nn5+nn4, sizeof(int));
958 for(i=0; i<kk5+kk4; i++)
959 J28[i] = calloc(kk5+kk4, sizeof(int));
960 //appendBlockMatrix(GG4, GG5, G32, nn4, nn4, nn5, nn5);
961 addSubMatrix(GG4, G28, kk4, nn4, 0, 0, 0, 0);
962 addSubMatrix(GG5, G28, kk5, nn5, 0, 0, kk4, nn4);
963 addSubMatrix(GG4, G28, nn4-kk4, nn4, kk4, 0, kk4+kk5, 0); // put the out-of-subspace vectors at the end
964 addSubMatrix(GG5, G28, nn5-kk5, nn5, kk5, 0, kk4+kk5+nn4-kk4, nn5);
965 //appendBlockMatrix(GGBar4, GGBar5, GBar32, nn4, nn4, nn5, nn5);
966 addSubMatrix(GGBar4, GBar28, kk4, nn4, 0, 0, 0, 0);
967 addSubMatrix(GGBar5, GBar28, kk5, nn5, 0, 0, kk4, nn4);
968 addSubMatrix(GGBar4, GBar28, nn4-kk4, nn4, kk4, 0, kk4+kk5, 0); // put the out-of-subspace vectors at the end
969 addSubMatrix(GGBar5, GBar28, nn5-kk5, nn5, kk5, 0, kk4+kk5+nn4-kk4, nn5);
970 appendVector(hh4, hh5, h28, nn4, nn5);
972 appendVector(DD4, DD5, D28, kk4, kk5);
973 appendBlockMatrix(JJ4, JJ5, J28, kk4, kk4, kk5, kk5);
976 int n29 = nn5+nn5; int k29 = kk5+kk5; int **G29; int **GBar29; int *h29; int Q29; int *D29; int **J29;
977 G29 = calloc(nn5+nn5,sizeof(int*));
978 GBar29 = calloc(nn5+nn5,sizeof(int*));
979 h29 = calloc(nn5+nn5,sizeof(int));
980 D29 = calloc(kk5+kk5,sizeof(int));
981 J29 = calloc(kk5+kk5,sizeof(int*));
982 for(i=0; i<nn5+nn5; i++) {
983 G29[i] = calloc(nn5+nn5, sizeof(int)); GBar29[i] = calloc(nn5+nn5, sizeof(int));
985 for(i=0; i<kk5+kk5; i++)
986 J29[i] = calloc(kk5+kk5, sizeof(int));
987 //appendBlockMatrix(GG5, GG5, G33, nn5, nn5, nn5, nn5);
988 addSubMatrix(GG5, G29, kk5, nn5, 0, 0, 0, 0);
989 addSubMatrix(GG5, G29, kk5, nn5, 0, 0, kk5, nn5);
990 addSubMatrix(GG5, G29, nn5-kk5, nn5, kk5, 0, kk5+kk5, 0); // put the out-of-subspace vectors at the end
991 addSubMatrix(GG5, G29, nn5-kk5, nn5, kk5, 0, kk5+kk5+nn5-kk5, nn5);
992 //appendBlockMatrix(GGBar5, GGBar5, GBar33, nn5, nn5, nn5, nn5);
993 addSubMatrix(GGBar5, GBar29, kk5, nn5, 0, 0, 0, 0);
994 addSubMatrix(GGBar5, GBar29, kk5, nn5, 0, 0, kk5, nn5);
995 addSubMatrix(GGBar5, GBar29, nn5-kk5, nn5, kk5, 0, kk5+kk5, 0); // put the out-of-subspace vectors at the end
996 addSubMatrix(GGBar5, GBar29, nn5-kk5, nn5, kk5, 0, kk5+kk5+nn5-kk5, nn5);
997 appendVector(hh5, hh5, h29, nn5, nn5);
999 appendVector(DD5, DD5, D29, kk5, kk5);
1000 appendBlockMatrix(JJ5, JJ5, J29, kk5, kk5, kk5, kk5);
1003 int n30 = nn5+nn6; int k30 = kk5+kk6; int **G30; int **GBar30; int *h30; int Q30; int *D30; int **J30;
1004 G30 = calloc(nn5+nn6,sizeof(int*));
1005 GBar30 = calloc(nn5+nn6,sizeof(int*));
1006 h30 = calloc(nn5+nn6,sizeof(int));
1007 D30 = calloc(kk5+kk6,sizeof(int));
1008 J30 = calloc(kk5+kk6,sizeof(int*));
1009 for(i=0; i<nn5+nn6; i++) {
1010 G30[i] = calloc(nn5+nn6, sizeof(int)); GBar30[i] = calloc(nn5+nn6, sizeof(int));
1012 for(i=0; i<kk5+kk6; i++)
1013 J30[i] = calloc(kk5+kk6, sizeof(int));
1014 //appendBlockMatrix(GG6, GG5, G34, nn6, nn6, nn5, nn5);
1015 addSubMatrix(GG6, G30, kk6, nn6, 0, 0, 0, 0);
1016 addSubMatrix(GG5, G30, kk5, nn5, 0, 0, kk6, nn6);
1017 addSubMatrix(GG6, G30, nn6-kk6, nn6, kk6, 0, kk6+kk5, 0); // put the out-of-subspace vectors at the end
1018 addSubMatrix(GG5, G30, nn5-kk5, nn5, kk5, 0, kk6+kk5+nn6-kk6, nn5);
1019 //appendBlockMatrix(GGBar6, GGBar5, GBar34, nn6, nn6, nn5, nn5);
1020 addSubMatrix(GGBar6, GBar30, kk6, nn6, 0, 0, 0, 0);
1021 addSubMatrix(GGBar5, GBar30, kk5, nn5, 0, 0, kk6, nn6);
1022 addSubMatrix(GGBar6, GBar30, nn6-kk6, nn6, kk6, 0, kk6+kk5, 0); // put the out-of-subspace vectors at the end
1023 addSubMatrix(GGBar5, GBar30, nn5-kk5, nn5, kk5, 0, kk6+kk5+nn6-kk6, nn5);
1024 appendVector(hh6, hh5, h30, nn6, nn5);
1026 appendVector(DD6, DD5, D30, kk6, kk5);
1027 appendBlockMatrix(JJ6, JJ5, J30, kk6, kk6, kk5, kk5);
1030 int n31 = nn5+nn7; int k31 = kk5+kk7; int **G31; int **GBar31; int *h31; int Q31; int *D31; int **J31;
1031 G31 = calloc(nn5+nn7,sizeof(int*));
1032 GBar31 = calloc(nn5+nn7,sizeof(int*));
1033 h31 = calloc(nn5+nn7,sizeof(int));
1034 D31 = calloc(kk5+kk7,sizeof(int));
1035 J31 = calloc(kk5+kk7,sizeof(int*));
1036 for(i=0; i<nn5+nn7; i++) {
1037 G31[i] = calloc(nn5+nn7, sizeof(int)); GBar31[i] = calloc(nn5+nn7, sizeof(int));
1039 for(i=0; i<kk5+kk7; i++)
1040 J31[i] = calloc(kk5+kk7, sizeof(int));
1041 //appendBlockMatrix(GG7, GG5, G35, nn7, nn7, nn5, nn5);
1042 addSubMatrix(GG7, G31, kk7, nn7, 0, 0, 0, 0);
1043 addSubMatrix(GG5, G31, kk5, nn5, 0, 0, kk7, nn7);
1044 addSubMatrix(GG7, G31, nn7-kk7, nn7, kk7, 0, kk7+kk5, 0); // put the out-of-subspace vectors at the end
1045 addSubMatrix(GG5, G31, nn5-kk5, nn5, kk5, 0, kk7+kk5+nn7-kk7, nn5);
1046 //appendBlockMatrix(GGBar7, GGBar5, GBar35, nn7, nn7, nn5, nn5);
1047 addSubMatrix(GGBar7, GBar31, kk7, nn7, 0, 0, 0, 0);
1048 addSubMatrix(GGBar5, GBar31, kk5, nn5, 0, 0, kk7, nn7);
1049 addSubMatrix(GGBar7, GBar31, nn7-kk7, nn7, kk7, 0, kk7+kk5, 0); // put the out-of-subspace vectors at the end
1050 addSubMatrix(GGBar5, GBar31, nn5-kk5, nn5, kk5, 0, kk7+kk5+nn7-kk7, nn5);
1051 appendVector(hh7, hh5, h31, nn7, nn5);
1053 appendVector(DD7, DD5, D31, kk7, kk5);
1054 appendBlockMatrix(JJ7, JJ5, J31, kk7, kk7, kk5, kk5);
1057 int n32 = nn6+nn1; int k32 = kk6+kk1; int **G32; int **GBar32; int *h32; int Q32; int *D32; int **J32;
1058 G32 = calloc(nn6+nn1,sizeof(int*));
1059 GBar32 = calloc(nn6+nn1,sizeof(int*));
1060 h32 = calloc(nn6+nn1,sizeof(int));
1061 D32 = calloc(kk6+kk1,sizeof(int));
1062 J32 = calloc(kk6+kk1,sizeof(int*));
1063 for(i=0; i<nn6+nn1; i++) {
1064 G32[i] = calloc(nn6+nn1, sizeof(int)); GBar32[i] = calloc(nn6+nn1, sizeof(int));
1066 for(i=0; i<kk6+kk1; i++)
1067 J32[i] = calloc(kk6+kk1, sizeof(int));
1068 //appendBlockMatrix(GG1, GG6, G36, nn1, nn1, nn6, nn6);
1069 addSubMatrix(GG1, G32, kk1, nn1, 0, 0, 0, 0);
1070 addSubMatrix(GG6, G32, kk6, nn6, 0, 0, kk1, nn1);
1071 addSubMatrix(GG1, G32, nn1-kk1, nn1, kk1, 0, kk1+kk6, 0); // put the out-of-subspace vectors at the end
1072 addSubMatrix(GG6, G32, nn6-kk6, nn6, kk6, 0, kk1+kk6+nn1-kk1, nn6);
1073 //appendBlockMatrix(GGBar1, GGBar6, GBar36, nn1, nn1, nn6, nn6);
1074 addSubMatrix(GGBar1, GBar32, kk1, nn1, 0, 0, 0, 0);
1075 addSubMatrix(GGBar6, GBar32, kk6, nn6, 0, 0, kk1, nn1);
1076 addSubMatrix(GGBar1, GBar32, nn1-kk1, nn1, kk1, 0, kk1+kk6, 0); // put the out-of-subspace vectors at the end
1077 addSubMatrix(GGBar6, GBar32, nn6-kk6, nn6, kk6, 0, kk1+kk6+nn1-kk1, nn6);
1078 appendVector(hh1, hh6, h32, nn1, nn6);
1080 appendVector(DD1, DD6, D32, kk1, kk6);
1081 appendBlockMatrix(JJ1, JJ6, J32, kk1, kk1, kk6, kk6);
1084 int n33 = nn6+nn2; int k33 = kk6+kk2; int **G33; int **GBar33; int *h33; int Q33; int *D33; int **J33;
1085 G33 = calloc(nn6+nn2,sizeof(int*));
1086 GBar33 = calloc(nn6+nn2,sizeof(int*));
1087 h33 = calloc(nn6+nn2,sizeof(int));
1088 D33 = calloc(kk6+kk2,sizeof(int));
1089 J33 = calloc(kk6+kk2,sizeof(int*));
1090 for(i=0; i<nn6+nn2; i++) {
1091 G33[i] = calloc(nn6+nn2, sizeof(int)); GBar33[i] = calloc(nn6+nn2, sizeof(int));
1093 for(i=0; i<kk6+kk2; i++)
1094 J33[i] = calloc(kk6+kk2, sizeof(int));
1095 //appendBlockMatrix(GG2, GG6, G37, nn2, nn2, nn6, nn6);
1096 addSubMatrix(GG2, G33, kk2, nn2, 0, 0, 0, 0);
1097 addSubMatrix(GG6, G33, kk6, nn6, 0, 0, kk2, nn2);
1098 addSubMatrix(GG2, G33, nn2-kk2, nn2, kk2, 0, kk2+kk6, 0); // put the out-of-subspace vectors at the end
1099 addSubMatrix(GG6, G33, nn6-kk6, nn6, kk6, 0, kk2+kk6+nn2-kk2, nn6);
1100 //appendBlockMatrix(GGBar2, GGBar6, GBar37, nn2, nn2, nn6, nn6);
1101 addSubMatrix(GGBar2, GBar33, kk2, nn2, 0, 0, 0, 0);
1102 addSubMatrix(GGBar6, GBar33, kk6, nn6, 0, 0, kk2, nn2);
1103 addSubMatrix(GGBar2, GBar33, nn2-kk2, nn2, kk2, 0, kk2+kk6, 0); // put the out-of-subspace vectors at the end
1104 addSubMatrix(GGBar6, GBar33, nn6-kk6, nn6, kk6, 0, kk2+kk6+nn2-kk2, nn6);
1105 appendVector(hh2, hh6, h33, nn2, nn6);
1107 appendVector(DD2, DD6, D33, kk2, kk6);
1108 appendBlockMatrix(JJ2, JJ6, J33, kk2, kk2, kk6, kk6);
1111 int n34 = nn6+nn3; int k34 = kk6+kk3; int **G34; int **GBar34; int *h34; int Q34; int *D34; int **J34;
1112 G34 = calloc(nn6+nn3,sizeof(int*));
1113 GBar34 = calloc(nn6+nn3,sizeof(int*));
1114 h34 = calloc(nn6+nn3,sizeof(int));
1115 D34 = calloc(kk6+kk3,sizeof(int));
1116 J34 = calloc(kk6+kk3,sizeof(int*));
1117 for(i=0; i<nn6+nn3; i++) {
1118 G34[i] = calloc(nn6+nn3, sizeof(int)); GBar34[i] = calloc(nn6+nn3, sizeof(int));
1120 for(i=0; i<kk6+kk3; i++)
1121 J34[i] = calloc(kk6+kk3, sizeof(int));
1122 //appendBlockMatrix(GG3, GG6, G38, nn3, nn3, nn6, nn6);
1123 addSubMatrix(GG3, G34, kk3, nn3, 0, 0, 0, 0);
1124 addSubMatrix(GG6, G34, kk6, nn6, 0, 0, kk3, nn3);
1125 addSubMatrix(GG3, G34, nn3-kk3, nn3, kk3, 0, kk3+kk6, 0); // put the out-of-subspace vectors at the end
1126 addSubMatrix(GG6, G34, nn6-kk6, nn6, kk6, 0, kk3+kk6+nn3-kk3, nn6);
1127 //appendBlockMatrix(GGBar3, GGBar6, GBar38, nn3, nn3, nn6, nn6);
1128 addSubMatrix(GGBar3, GBar34, kk3, nn3, 0, 0, 0, 0);
1129 addSubMatrix(GGBar6, GBar34, kk6, nn6, 0, 0, kk3, nn3);
1130 addSubMatrix(GGBar3, GBar34, nn3-kk3, nn3, kk3, 0, kk3+kk6, 0); // put the out-of-subspace vectors at the end
1131 addSubMatrix(GGBar6, GBar34, nn6-kk6, nn6, kk6, 0, kk3+kk6+nn3-kk3, nn6);
1132 appendVector(hh3, hh6, h34, nn3, nn6);
1134 appendVector(DD3, DD6, D34, kk3, kk6);
1135 appendBlockMatrix(JJ3, JJ6, J34, kk3, kk3, kk6, kk6);
1138 int n35 = nn6+nn4; int k35 = kk6+kk4; int **G35; int **GBar35; int *h35; int Q35; int *D35; int **J35;
1139 G35 = calloc(nn6+nn4,sizeof(int*));
1140 GBar35 = calloc(nn6+nn4,sizeof(int*));
1141 h35 = calloc(nn6+nn4,sizeof(int));
1142 D35 = calloc(kk6+kk4,sizeof(int));
1143 J35 = calloc(kk6+kk4,sizeof(int*));
1144 for(i=0; i<nn6+nn4; i++) {
1145 G35[i] = calloc(nn6+nn4, sizeof(int)); GBar35[i] = calloc(nn6+nn4, sizeof(int));
1147 for(i=0; i<kk6+kk4; i++)
1148 J35[i] = calloc(kk6+kk4, sizeof(int));
1149 //appendBlockMatrix(GG4, GG6, G39, nn4, nn4, nn6, nn6);
1150 addSubMatrix(GG4, G35, kk4, nn4, 0, 0, 0, 0);
1151 addSubMatrix(GG6, G35, kk6, nn6, 0, 0, kk4, nn4);
1152 addSubMatrix(GG4, G35, nn4-kk4, nn4, kk4, 0, kk4+kk6, 0); // put the out-of-subspace vectors at the end
1153 addSubMatrix(GG6, G35, nn6-kk6, nn6, kk6, 0, kk4+kk6+nn4-kk4, nn6);
1154 //appendBlockMatrix(GGBar4, GGBar6, GBar39, nn4, nn4, nn6, nn6);
1155 addSubMatrix(GGBar4, GBar35, kk4, nn4, 0, 0, 0, 0);
1156 addSubMatrix(GGBar6, GBar35, kk6, nn6, 0, 0, kk4, nn4);
1157 addSubMatrix(GGBar4, GBar35, nn4-kk4, nn4, kk4, 0, kk4+kk6, 0); // put the out-of-subspace vectors at the end
1158 addSubMatrix(GGBar6, GBar35, nn6-kk6, nn6, kk6, 0, kk4+kk6+nn4-kk4, nn6);
1159 appendVector(hh4, hh6, h35, nn4, nn6);
1161 appendVector(DD4, DD6, D35, kk4, kk6);
1162 appendBlockMatrix(JJ4, JJ6, J35, kk4, kk4, kk6, kk6);
1165 int n36 = nn6+nn5; int k36 = kk6+kk5; int **G36; int **GBar36; int *h36; int Q36; int *D36; int **J36;
1166 G36 = calloc(nn6+nn5,sizeof(int*));
1167 GBar36 = calloc(nn6+nn5,sizeof(int*));
1168 h36 = calloc(nn6+nn5,sizeof(int));
1169 D36 = calloc(kk6+kk5,sizeof(int));
1170 J36 = calloc(kk6+kk5,sizeof(int*));
1171 for(i=0; i<nn6+nn5; i++) {
1172 G36[i] = calloc(nn6+nn5, sizeof(int)); GBar36[i] = calloc(nn6+nn5, sizeof(int));
1174 for(i=0; i<kk6+kk5; i++)
1175 J36[i] = calloc(kk6+kk5, sizeof(int));
1176 //appendBlockMatrix(GG5, GG6, G40, nn5, nn5, nn6, nn6);
1177 addSubMatrix(GG5, G36, kk5, nn5, 0, 0, 0, 0);
1178 addSubMatrix(GG6, G36, kk6, nn6, 0, 0, kk5, nn5);
1179 addSubMatrix(GG5, G36, nn5-kk5, nn5, kk5, 0, kk5+kk6, 0); // put the out-of-subspace vectors at the end
1180 addSubMatrix(GG6, G36, nn6-kk6, nn6, kk6, 0, kk5+kk6+nn5-kk5, nn6);
1181 //appendBlockMatrix(GGBar5, GGBar6, GBar40, nn5, nn5, nn6, nn6);
1182 addSubMatrix(GGBar5, GBar36, kk5, nn5, 0, 0, 0, 0);
1183 addSubMatrix(GGBar6, GBar36, kk6, nn6, 0, 0, kk5, nn5);
1184 addSubMatrix(GGBar5, GBar36, nn5-kk5, nn5, kk5, 0, kk5+kk6, 0); // put the out-of-subspace vectors at the end
1185 addSubMatrix(GGBar6, GBar36, nn6-kk6, nn6, kk6, 0, kk5+kk6+nn5-kk5, nn6);
1186 appendVector(hh5, hh6, h36, nn5, nn6);
1188 appendVector(DD5, DD6, D36, kk5, kk6);
1189 appendBlockMatrix(JJ5, JJ6, J36, kk5, kk5, kk6, kk6);
1191 // phiprime * phiprime
1192 int n37 = nn6+nn6; int k37 = kk6+kk6; int **G37; int **GBar37; int *h37; int Q37; int *D37; int **J37;
1193 G37 = calloc(nn6+nn6,sizeof(int*));
1194 GBar37 = calloc(nn6+nn6,sizeof(int*));
1195 h37 = calloc(nn6+nn6,sizeof(int));
1196 D37 = calloc(kk6+kk6,sizeof(int));
1197 J37 = calloc(kk6+kk6,sizeof(int*));
1198 for(i=0; i<nn6+nn6; i++) {
1199 G37[i] = calloc(nn6+nn6, sizeof(int)); GBar37[i] = calloc(nn6+nn6, sizeof(int));
1201 for(i=0; i<kk6+kk6; i++)
1202 J37[i] = calloc(kk6+kk6, sizeof(int));
1203 //appendBlockMatrix(GG6, GG6, G41, nn6, nn6, nn6, nn6);
1204 addSubMatrix(GG6, G37, kk6, nn6, 0, 0, 0, 0);
1205 addSubMatrix(GG6, G37, kk6, nn6, 0, 0, kk6, nn6);
1206 addSubMatrix(GG6, G37, nn6-kk6, nn6, kk6, 0, kk6+kk6, 0); // put the out-of-subspace vectors at the end
1207 addSubMatrix(GG6, G37, nn6-kk6, nn6, kk6, 0, kk6+kk6+nn6-kk6, nn6);
1208 //appendBlockMatrix(GGBar6, GGBar6, GBar41, nn6, nn6, nn6, nn6);
1209 addSubMatrix(GGBar6, GBar37, kk6, nn6, 0, 0, 0, 0);
1210 addSubMatrix(GGBar6, GBar37, kk6, nn6, 0, 0, kk6, nn6);
1211 addSubMatrix(GGBar6, GBar37, nn6-kk6, nn6, kk6, 0, kk6+kk6, 0); // put the out-of-subspace vectors at the end
1212 addSubMatrix(GGBar6, GBar37, nn6-kk6, nn6, kk6, 0, kk6+kk6+nn6-kk6, nn6);
1213 appendVector(hh6, hh6, h37, nn6, nn6);
1215 appendVector(DD6, DD6, D37, kk6, kk6);
1216 appendBlockMatrix(JJ6, JJ6, J37, kk6, kk6, kk6, kk6);
1218 // phidprime * phiprime
1219 int n38 = nn6+nn7; int k38 = kk6+kk7; int **G38; int **GBar38; int *h38; int Q38; int *D38; int **J38;
1220 G38 = calloc(nn6+nn7,sizeof(int*));
1221 GBar38 = calloc(nn6+nn7,sizeof(int*));
1222 h38 = calloc(nn6+nn7,sizeof(int));
1223 D38 = calloc(kk6+kk7,sizeof(int));
1224 J38 = calloc(kk6+kk7,sizeof(int*));
1225 for(i=0; i<nn6+nn7; i++) {
1226 G38[i] = calloc(nn6+nn7, sizeof(int)); GBar38[i] = calloc(nn6+nn7, sizeof(int));
1228 for(i=0; i<kk6+kk7; i++)
1229 J38[i] = calloc(kk6+kk7, sizeof(int));
1230 //appendBlockMatrix(GG7, GG6, G42, nn7, nn7, nn6, nn6);
1231 addSubMatrix(GG7, G38, kk7, nn7, 0, 0, 0, 0);
1232 addSubMatrix(GG6, G38, kk6, nn6, 0, 0, kk7, nn7);
1233 addSubMatrix(GG7, G38, nn7-kk7, nn7, kk7, 0, kk7+kk6, 0); // put the out-of-subspace vectors at the end
1234 addSubMatrix(GG6, G38, nn6-kk6, nn6, kk6, 0, kk7+kk6+nn7-kk7, nn6);
1235 //appendBlockMatrix(GGBar7, GGBar6, GBar42, nn7, nn7, nn6, nn6);
1236 addSubMatrix(GGBar7, GBar38, kk7, nn7, 0, 0, 0, 0);
1237 addSubMatrix(GGBar6, GBar38, kk6, nn6, 0, 0, kk7, nn7);
1238 addSubMatrix(GGBar7, GBar38, nn7-kk7, nn7, kk7, 0, kk7+kk6, 0); // put the out-of-subspace vectors at the end
1239 addSubMatrix(GGBar6, GBar38, nn6-kk6, nn6, kk6, 0, kk7+kk6+nn7-kk7, nn6);
1240 appendVector(hh7, hh6, h38, nn7, nn6);
1242 appendVector(DD7, DD6, D38, kk7, kk6);
1243 appendBlockMatrix(JJ7, JJ6, J38, kk7, kk7, kk6, kk6);
1246 int n39 = nn7+nn1; int k39 = kk7+kk1; int **G39; int **GBar39; int *h39; int Q39; int *D39; int **J39;
1247 G39 = calloc(nn7+nn1,sizeof(int*));
1248 GBar39 = calloc(nn7+nn1,sizeof(int*));
1249 h39 = calloc(nn7+nn1,sizeof(int));
1250 D39 = calloc(kk7+kk1,sizeof(int));
1251 J39 = calloc(kk7+kk1,sizeof(int*));
1252 for(i=0; i<nn7+nn1; i++) {
1253 G39[i] = calloc(nn7+nn1, sizeof(int)); GBar39[i] = calloc(nn7+nn1, sizeof(int));
1255 for(i=0; i<kk7+kk1; i++)
1256 J39[i] = calloc(kk7+kk1, sizeof(int));
1257 //appendBlockMatrix(GG1, GG7, G43, nn1, nn1, nn7, nn7);
1258 addSubMatrix(GG1, G39, kk1, nn1, 0, 0, 0, 0);
1259 addSubMatrix(GG7, G39, kk7, nn7, 0, 0, kk1, nn1);
1260 addSubMatrix(GG1, G39, nn1-kk1, nn1, kk1, 0, kk1+kk7, 0); // put the out-of-subspace vectors at the end
1261 addSubMatrix(GG7, G39, nn7-kk7, nn7, kk7, 0, kk1+kk7+nn1-kk1, nn7);
1262 //appendBlockMatrix(GGBar1, GGBar7, GBar43, nn1, nn1, nn7, nn7);
1263 addSubMatrix(GGBar1, GBar39, kk1, nn1, 0, 0, 0, 0);
1264 addSubMatrix(GGBar7, GBar39, kk7, nn7, 0, 0, kk1, nn1);
1265 addSubMatrix(GGBar1, GBar39, nn1-kk1, nn1, kk1, 0, kk1+kk7, 0); // put the out-of-subspace vectors at the end
1266 addSubMatrix(GGBar7, GBar39, nn7-kk7, nn7, kk7, 0, kk1+kk7+nn1-kk1, nn7);
1267 appendVector(hh1, hh7, h39, nn1, nn7);
1269 appendVector(DD1, DD7, D39, kk1, kk7);
1270 appendBlockMatrix(JJ1, JJ7, J39, kk1, kk1, kk7, kk7);
1273 int n40 = nn7+nn2; int k40 = kk7+kk2; int **G40; int **GBar40; int *h40; int Q40; int *D40; int **J40;
1274 G40 = calloc(nn7+nn2,sizeof(int*));
1275 GBar40 = calloc(nn7+nn2,sizeof(int*));
1276 h40 = calloc(nn7+nn2,sizeof(int));
1277 D40 = calloc(kk7+kk2,sizeof(int));
1278 J40 = calloc(kk7+kk2,sizeof(int*));
1279 for(i=0; i<nn7+nn2; i++) {
1280 G40[i] = calloc(nn7+nn2, sizeof(int)); GBar40[i] = calloc(nn7+nn2, sizeof(int));
1282 for(i=0; i<kk7+kk2; i++)
1283 J40[i] = calloc(kk7+kk2, sizeof(int));
1284 //appendBlockMatrix(GG2, GG7, G44, nn2, nn2, nn7, nn7);
1285 addSubMatrix(GG2, G40, kk2, nn2, 0, 0, 0, 0);
1286 addSubMatrix(GG7, G40, kk7, nn7, 0, 0, kk2, nn2);
1287 addSubMatrix(GG1, G40, nn2-kk2, nn2, kk2, 0, kk2+kk7, 0); // put the out-of-subspace vectors at the end
1288 addSubMatrix(GG7, G40, nn7-kk7, nn7, kk7, 0, kk2+kk7+nn2-kk2, nn7);
1289 //appendBlockMatrix(GGBar2, GGBar7, GBar44, nn2, nn2, nn7, nn7);
1290 addSubMatrix(GGBar2, GBar40, kk2, nn2, 0, 0, 0, 0);
1291 addSubMatrix(GGBar7, GBar40, kk7, nn7, 0, 0, kk2, nn2);
1292 addSubMatrix(GGBar1, GBar40, nn2-kk2, nn2, kk2, 0, kk2+kk7, 0); // put the out-of-subspace vectors at the end
1293 addSubMatrix(GGBar7, GBar40, nn7-kk7, nn7, kk7, 0, kk2+kk7+nn2-kk2, nn7);
1294 appendVector(hh2, hh7, h40, nn2, nn7);
1296 appendVector(DD2, DD7, D40, kk2, kk7);
1297 appendBlockMatrix(JJ2, JJ7, J40, kk2, kk2, kk7, kk7);
1300 int n41 = nn7+nn3; int k41 = kk7+kk3; int **G41; int **GBar41; int *h41; int Q41; int *D41; int **J41;
1301 G41 = calloc(nn7+nn3,sizeof(int*));
1302 GBar41 = calloc(nn7+nn3,sizeof(int*));
1303 h41 = calloc(nn7+nn3,sizeof(int));
1304 D41 = calloc(kk7+kk3,sizeof(int));
1305 J41 = calloc(kk7+kk3,sizeof(int*));
1306 for(i=0; i<nn7+nn3; i++) {
1307 G41[i] = calloc(nn7+nn3, sizeof(int)); GBar41[i] = calloc(nn7+nn3, sizeof(int));
1309 for(i=0; i<kk7+kk3; i++)
1310 J41[i] = calloc(kk7+kk3, sizeof(int));
1311 //appendBlockMatrix(GG3, GG7, G45, nn3, nn3, nn7, nn7);
1312 addSubMatrix(GG3, G41, kk3, nn3, 0, 0, 0, 0);
1313 addSubMatrix(GG7, G41, kk7, nn7, 0, 0, kk3, nn3);
1314 addSubMatrix(GG3, G41, nn3-kk3, nn3, kk3, 0, kk3+kk7, 0); // put the out-of-subspace vectors at the end
1315 addSubMatrix(GG7, G41, nn7-kk7, nn7, kk7, 0, kk3+kk7+nn3-kk3, nn7);
1316 //appendBlockMatrix(GGBar3, GGBar7, GBar45, nn3, nn3, nn7, nn7);
1317 addSubMatrix(GGBar3, GBar41, kk3, nn3, 0, 0, 0, 0);
1318 addSubMatrix(GGBar7, GBar41, kk7, nn7, 0, 0, kk3, nn3);
1319 addSubMatrix(GGBar3, GBar41, nn3-kk3, nn3, kk3, 0, kk3+kk7, 0); // put the out-of-subspace vectors at the end
1320 addSubMatrix(GGBar7, GBar41, nn7-kk7, nn7, kk7, 0, kk3+kk7+nn3-kk3, nn7);
1321 appendVector(hh3, hh7, h41, nn3, nn7);
1323 appendVector(DD3, DD7, D41, kk3, kk7);
1324 appendBlockMatrix(JJ3, JJ7, J41, kk3, kk3, kk7, kk7);
1327 int n42 = nn7+nn4; int k42 = kk7+kk4; int **G42; int **GBar42; int *h42; int Q42; int *D42; int **J42;
1328 G42 = calloc(nn7+nn4,sizeof(int*));
1329 GBar42 = calloc(nn7+nn4,sizeof(int*));
1330 h42 = calloc(nn7+nn4,sizeof(int));
1331 D42 = calloc(kk7+kk4,sizeof(int));
1332 J42 = calloc(kk7+kk4,sizeof(int*));
1333 for(i=0; i<nn7+nn4; i++) {
1334 G42[i] = calloc(nn7+nn4, sizeof(int)); GBar42[i] = calloc(nn7+nn4, sizeof(int));
1336 for(i=0; i<kk7+kk4; i++)
1337 J42[i] = calloc(kk7+kk4, sizeof(int));
1338 //appendBlockMatrix(GG4, GG7, G46, nn4, nn4, nn7, nn7);
1339 addSubMatrix(GG4, G42, kk4, nn4, 0, 0, 0, 0);
1340 addSubMatrix(GG7, G42, kk7, nn7, 0, 0, kk4, nn4);
1341 addSubMatrix(GG4, G42, nn4-kk4, nn4, kk4, 0, kk4+kk7, 0); // put the out-of-subspace vectors at the end
1342 addSubMatrix(GG7, G42, nn7-kk7, nn7, kk7, 0, kk4+kk7+nn4-kk4, nn7);
1343 //appendBlockMatrix(GGBar4, GGBar7, GBar46, nn4, nn4, nn7, nn7);
1344 addSubMatrix(GGBar4, GBar42, kk4, nn4, 0, 0, 0, 0);
1345 addSubMatrix(GGBar7, GBar42, kk7, nn7, 0, 0, kk4, nn4);
1346 addSubMatrix(GGBar4, GBar42, nn4-kk4, nn4, kk4, 0, kk4+kk7, 0); // put the out-of-subspace vectors at the end
1347 addSubMatrix(GGBar7, GBar42, nn7-kk7, nn7, kk7, 0, kk4+kk7+nn4-kk4, nn7);
1348 appendVector(hh4, hh7, h42, nn4, nn7);
1350 appendVector(DD4, DD7, D42, kk4, kk7);
1351 appendBlockMatrix(JJ4, JJ7, J42, kk4, kk4, kk7, kk7);
1354 int n43 = nn7+nn5; int k43 = kk7+kk5; int **G43; int **GBar43; int *h43; int Q43; int *D43; int **J43;
1355 G43 = calloc(nn7+nn5,sizeof(int*));
1356 GBar43 = calloc(nn7+nn5,sizeof(int*));
1357 h43 = calloc(nn7+nn5,sizeof(int));
1358 D43 = calloc(kk7+kk5,sizeof(int));
1359 J43 = calloc(kk7+kk5,sizeof(int*));
1360 for(i=0; i<nn7+nn5; i++) {
1361 G43[i] = calloc(nn7+nn5, sizeof(int)); GBar43[i] = calloc(nn7+nn5, sizeof(int));
1363 for(i=0; i<kk7+kk5; i++)
1364 J43[i] = calloc(kk7+kk5, sizeof(int));
1365 //appendBlockMatrix(GG5, GG7, G47, nn5, nn5, nn7, nn7);
1366 addSubMatrix(GG5, G43, kk5, nn5, 0, 0, 0, 0);
1367 addSubMatrix(GG7, G43, kk7, nn7, 0, 0, kk5, nn5);
1368 addSubMatrix(GG5, G43, nn5-kk5, nn5, kk5, 0, kk5+kk7, 0); // put the out-of-subspace vectors at the end
1369 addSubMatrix(GG7, G43, nn7-kk7, nn7, kk7, 0, kk5+kk7+nn5-kk5, nn7);
1370 //appendBlockMatrix(GGBar5, GGBar7, GBar47, nn5, nn5, nn7, nn7);
1371 addSubMatrix(GGBar5, GBar43, kk5, nn5, 0, 0, 0, 0);
1372 addSubMatrix(GGBar7, GBar43, kk7, nn7, 0, 0, kk5, nn5);
1373 addSubMatrix(GGBar5, GBar43, nn5-kk5, nn5, kk5, 0, kk5+kk7, 0); // put the out-of-subspace vectors at the end
1374 addSubMatrix(GGBar7, GBar43, nn7-kk7, nn7, kk7, 0, kk5+kk7+nn5-kk5, nn7);
1375 appendVector(hh5, hh7, h43, nn5, nn7);
1377 appendVector(DD5, DD7, D43, kk5, kk7);
1378 appendBlockMatrix(JJ5, JJ7, J43, kk5, kk5, kk7, kk7);
1380 // phiprime * phidprime
1381 int n44 = nn7+nn6; int k44 = kk7+kk6; int **G44; int **GBar44; int *h44; int Q44; int *D44; int **J44;
1382 G44 = calloc(nn7+nn6,sizeof(int*));
1383 GBar44 = calloc(nn7+nn6,sizeof(int*));
1384 h44 = calloc(nn7+nn6,sizeof(int));
1385 D44 = calloc(kk7+kk6,sizeof(int));
1386 J44 = calloc(kk7+kk6,sizeof(int*));
1387 for(i=0; i<nn7+nn6; i++) {
1388 G44[i] = calloc(nn7+nn6, sizeof(int)); GBar44[i] = calloc(nn7+nn6, sizeof(int));
1390 for(i=0; i<kk7+kk6; i++)
1391 J44[i] = calloc(kk7+kk6, sizeof(int));
1392 //appendBlockMatrix(GG6, GG7, G48, nn6, nn6, nn7, nn7);
1393 addSubMatrix(GG6, G44, kk6, nn6, 0, 0, 0, 0);
1394 addSubMatrix(GG7, G44, kk7, nn7, 0, 0, kk6, nn6);
1395 addSubMatrix(GG6, G44, nn6-kk6, nn6, kk6, 0, kk6+kk7, 0); // put the out-of-subspace vectors at the end
1396 addSubMatrix(GG7, G44, nn7-kk7, nn7, kk7, 0, kk6+kk7+nn6-kk6, nn7);
1397 //appendBlockMatrix(GGBar6, GGBar7, GBar48, nn6, nn6, nn7, nn7);
1398 addSubMatrix(GGBar6, GBar44, kk6, nn6, 0, 0, 0, 0);
1399 addSubMatrix(GGBar7, GBar44, kk7, nn7, 0, 0, kk6, nn6);
1400 addSubMatrix(GGBar6, GBar44, nn6-kk6, nn6, kk6, 0, kk6+kk7, 0); // put the out-of-subspace vectors at the end
1401 addSubMatrix(GGBar7, GBar44, nn7-kk7, nn7, kk7, 0, kk6+kk7+nn6-kk6, nn7);
1402 appendVector(hh6, hh7, h44, nn6, nn7);
1404 appendVector(DD6, DD7, D44, kk6, kk7);
1405 appendBlockMatrix(JJ6, JJ7, J44, kk6, kk6, kk7, kk7);
1407 // phidprime * phidprime
1408 int n45 = nn7+nn7; int k45 = kk7+kk7; int **G45; int **GBar45; int *h45; int Q45; int *D45; int **J45;
1409 G45 = calloc(nn7+nn7,sizeof(int*));
1410 GBar45 = calloc(nn7+nn7,sizeof(int*));
1411 h45 = calloc(nn7+nn7,sizeof(int));
1412 D45 = calloc(kk7+kk7,sizeof(int));
1413 J45 = calloc(kk7+kk7,sizeof(int*));
1414 for(i=0; i<nn7+nn7; i++) {
1415 G45[i] = calloc(nn7+nn7, sizeof(int)); GBar45[i] = calloc(nn7+nn7, sizeof(int));
1417 for(i=0; i<kk7+kk7; i++)
1418 J45[i] = calloc(kk7+kk7, sizeof(int));
1419 //appendBlockMatrix(GG7, GG7, G49, nn7, nn7, nn7, nn7);
1420 addSubMatrix(GG7, G45, kk7, nn7, 0, 0, 0, 0);
1421 addSubMatrix(GG7, G45, kk7, nn7, 0, 0, kk7, nn7);
1422 addSubMatrix(GG7, G45, nn7-kk7, nn7, kk7, 0, kk7+kk7, 0); // put the out-of-subspace vectors at the end
1423 addSubMatrix(GG7, G45, nn7-kk7, nn7, kk7, 0, kk7+kk7+nn7-kk7, nn7);
1424 //appendBlockMatrix(GGBar7, GGBar7, GBar49, nn7, nn7, nn7, nn7);
1425 addSubMatrix(GGBar7, GBar45, kk7, nn7, 0, 0, 0, 0);
1426 addSubMatrix(GGBar7, GBar45, kk7, nn7, 0, 0, kk7, nn7);
1427 addSubMatrix(GGBar7, GBar45, nn7-kk7, nn7, kk7, 0, kk7+kk7, 0); // put the out-of-subspace vectors at the end
1428 addSubMatrix(GGBar7, GBar45, nn7-kk7, nn7, kk7, 0, kk7+kk7+nn7-kk7, nn7);
1429 appendVector(hh7, hh7, h45, nn7, nn7);
1431 appendVector(DD7, DD7, D45, kk7, kk7);
1432 appendBlockMatrix(JJ7, JJ7, J45, kk7, kk7, kk7, kk7);
1434 // b60 * b66 + b66 * b60
1435 int n46 = 12; int k46 = 11;
1436 int (*(G46[])) = { (int[]) {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1437 int (*(GBar46[])) = { (int[]) {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};
1438 int h46[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
1440 int D46[] = {0, 0, 0, 0, 0, 4, 4, 4, 4, 4, 4};
1441 int (*(J46[])) = { (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1443 // e6 * o6 + e6 * o6
1444 int n47 = 12; int k47 = 11;
1445 int (*(G47[])) = { (int[]) {1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, (int[]) {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}};
1446 int (*(GBar47[])) = { (int[]) {0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}, (int[]) {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, (int[]) {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};
1447 int h47[] = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
1449 int D47[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
1450 int (*(J47[])) = { (int[]) {0, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}, (int[]) {4, 0, 4, 4, 4, 4, 4, 4, 4, 4, 4}, (int[]) {4, 4, 0, 4, 4, 4, 4, 4, 4, 4, 4}, (int[]) {4, 4, 4, 0, 4, 4, 4, 4, 4, 4, 4}, (int[]) {4, 4, 4, 4, 0, 4, 4, 4, 4, 4, 4}, (int[]) {4, 4, 4, 4, 4, 0, 4, 4, 4, 4, 4}, (int[]) {4, 4, 4, 4, 4, 4, 0, 4, 4, 4, 4}, (int[]) {4, 4, 4, 4, 4, 4, 4, 0, 4, 4, 4}, (int[]) {4, 4, 4, 4, 4, 4, 4, 4, 0, 4, 4}, (int[]) {4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 4}, (int[]) {4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0}};
1455 int *K; int ***G; int ***GBar; int **h; int *Q; int **D; int ***J;
1456 double complex Gamma[(int)pow(47,N/12)]; // prefactor in front of resultant state
1457 G = calloc(pow(47,N/12),sizeof(int*)); GBar = calloc(pow(47,N/12),sizeof(int*));
1458 h = calloc(pow(47,N/12),sizeof(int*));
1460 J = calloc(pow(47,N/12),sizeof(int*)); D = calloc(pow(47,N/12),sizeof(int*)); Q = calloc(pow(47,N/12),sizeof(int));
1462 K = calloc(pow(47,N/12), sizeof(int));
1464 int origK, origQ, *origD;
1466 int **origG, **origGBar;
1468 double complex origGamma;
1470 int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
1473 for(j=0; j<pow(47,N/12); j++) { // there will be 47^(N/12) combinations when using k=12 tensor factors
1479 for(k=0; k<N/12; k++) {
1480 K[j] += (((combination%47)==46)*k47 + ((combination%47)==45)*k46 + ((combination%47)==44)*k45 + ((combination%47)==43)*k44 + ((combination%47)==42)*k43 + ((combination%47)==41)*k42 + ((combination%47)==40)*k41 + ((combination%47)==39)*k40 +
1481 ((combination%47)==38)*k39 + ((combination%47)==37)*k38 + ((combination%47)==36)*k37 + ((combination%47)==35)*k36 + ((combination%47)==34)*k35 + ((combination%47)==33)*k34 + ((combination%47)==32)*k33 + ((combination%47)==31)*k32 + ((combination%47)==30)*k31 + ((combination%47)==29)*k30 +
1482 ((combination%47)==28)*k29 + ((combination%47)==27)*k28 + ((combination%47)==26)*k27 + ((combination%47)==25)*k26 + ((combination%47)==24)*k25 + ((combination%47)==23)*k24 + ((combination%47)==22)*k23 + ((combination%47)==21)*k22 + ((combination%47)==20)*k21 + ((combination%47)==19)*k20 +
1483 ((combination%47)==18)*k19 + ((combination%47)==17)*k18 + ((combination%47)==16)*k17 + ((combination%47)==15)*k16 + ((combination%47)==14)*k15 + ((combination%47)==13)*k14 + ((combination%47)==12)*k13 + ((combination%47)==11)*k12 + ((combination%47)==10)*k11 + ((combination%47)==9)*k10 +
1484 ((combination%47)==8)*k9 + ((combination%47)==7)*k8 + ((combination%47)==6)*k7 + ((combination%47)==5)*k6 + ((combination%47)==4)*k5 + ((combination%47)==3)*k4 + ((combination%47)==2)*k3 + ((combination%47)==1)*k2 + ((combination%47)==0)*k1);
1491 G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
1492 h[j] = calloc(N, sizeof(int));
1495 J[j] = calloc(K[j], sizeof(int*)); D[j] = calloc(K[j], sizeof(int));
1496 for(k=0; k<K[j]; k++)
1497 J[j][k] = calloc(K[j], sizeof(int));
1500 for(k=0; k<N; k++) {
1501 G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
1504 int Kcounter = 0; // Kcounter keeps track of the K<=N that we have added already to the G rows etc. for each combination that is indexed by the digits (base 3) of 'j' in that we go through with 'k'
1505 int Kcombo; // Kcombo stores the k<(n1=n2=n3) dimension of the member of the combination that we are currently adding
1506 for(k=0; k<N/12; k++) {
1508 Q[j] += (((combination%47)==46)*Q47 + ((combination%47)==45)*Q46 + ((combination%47)==44)*Q45 + ((combination%47)==43)*Q44 + ((combination%47)==42)*Q43 + ((combination%47)==41)*Q42 + ((combination%47)==40)*Q41 + ((combination%47)==39)*Q40 +
1509 ((combination%47)==38)*Q39 + ((combination%47)==37)*Q38 + ((combination%47)==36)*Q37 + ((combination%47)==35)*Q36 + ((combination%47)==34)*Q35 + ((combination%47)==33)*Q34 + ((combination%47)==32)*Q33 + ((combination%47)==31)*Q32 + ((combination%47)==30)*Q31 + ((combination%47)==29)*Q30 +
1510 ((combination%47)==28)*Q29 + ((combination%47)==27)*Q28 + ((combination%47)==26)*Q27 + ((combination%47)==25)*Q26 + ((combination%47)==24)*Q25 + ((combination%47)==23)*Q24 + ((combination%47)==22)*Q23 + ((combination%47)==21)*Q22 + ((combination%47)==20)*Q21 + ((combination%47)==19)*Q20 +
1511 ((combination%47)==18)*Q19 + ((combination%47)==17)*Q18 + ((combination%47)==16)*Q17 + ((combination%47)==15)*Q16 + ((combination%47)==14)*Q15 + ((combination%47)==13)*Q14 + ((combination%47)==12)*Q13 + ((combination%47)==11)*Q12 + ((combination%47)==10)*Q11 + ((combination%47)==9)*Q10 +
1512 ((combination%47)==8)*Q9 + ((combination%47)==7)*Q8 + ((combination%47)==6)*Q7 + ((combination%47)==5)*Q6 + ((combination%47)==4)*Q5 + ((combination%47)==3)*Q4 + ((combination%47)==2)*Q3 + ((combination%47)==1)*Q2 + ((combination%47)==0)*Q1);
1515 Gamma[j] *= (((combination%47)==46)*coeffe6*coeffo6*sqrt(2.0) + ((combination%47)==45)*coeffb60*coeffb66*sqrt(2.0) + ((combination%47)==44)*coeffphidprime*coeffphidprime + ((combination%47)==43)*coeffphiprime*coeffphidprime + ((combination%47)==42)*coeffk6*coeffphidprime + ((combination%47)==41)*coeffo6*coeffphidprime + ((combination%47)==40)*coeffe6*coeffphidprime + ((combination%47)==39)*coeffb66*coeffphidprime + ((combination%47)==38)*coeffb60*coeffphidprime
1516 + ((combination%47)==37)*coeffphidprime*coeffphiprime + ((combination%47)==36)*coeffphiprime*coeffphiprime + ((combination%47)==35)*coeffk6*coeffphiprime + ((combination%47)==34)*coeffo6*coeffphiprime + ((combination%47)==33)*coeffe6*coeffphiprime + ((combination%47)==32)*coeffb66*coeffphiprime + ((combination%47)==31)*coeffb60*coeffphiprime
1517 + ((combination%47)==30)*coeffphidprime*coeffk6 + ((combination%47)==29)*coeffphiprime*coeffk6 + ((combination%47)==28)*coeffk6*coeffk6 + ((combination%47)==27)*coeffo6*coeffk6 + ((combination%47)==26)*coeffe6*coeffk6 + ((combination%47)==25)*coeffb66*coeffk6 + ((combination%47)==24)*coeffb60*coeffk6
1518 + ((combination%47)==23)*coeffphidprime*coeffo6 + ((combination%47)==22)*coeffphiprime*coeffo6 + ((combination%47)==21)*coeffk6*coeffo6 + ((combination%47)==20)*coeffo6*coeffo6 + ((combination%47)==19)*coeffb66*coeffo6 + ((combination%47)==18)*coeffb60*coeffo6
1519 + ((combination%47)==17)*coeffphidprime*coeffe6 + ((combination%47)==16)*coeffphiprime*coeffe6 + ((combination%47)==15)*coeffk6*coeffe6 + ((combination%47)==14)*coeffe6*coeffe6 + ((combination%47)==13)*coeffb66*coeffe6 + ((combination%47)==12)*coeffb60*coeffe6
1520 + ((combination%47)==11)*coeffphidprime*coeffb66 + ((combination%47)==10)*coeffphiprime*coeffb66 + ((combination%47)==9)*coeffk6*coeffb66 + ((combination%47)==8)*coeffo6*coeffb66 + ((combination%47)==7)*coeffe6*coeffb66 + ((combination%47)==6)*coeffb66*coeffb66
1521 + ((combination%47)==5)*coeffphidprime*coeffb60 + ((combination%47)==4)*coeffphiprime*coeffb60 + ((combination%47)==3)*coeffk6*coeffb60 + ((combination%47)==2)*coeffo6*coeffb60 + ((combination%47)==1)*coeffe6*coeffb60 + ((combination%47)==0)*coeffb60*coeffb60);
1523 Kcombo = (((combination%47)==46)*k47 + ((combination%47)==45)*k46 + ((combination%47)==44)*k45 + ((combination%47)==43)*k44 + ((combination%47)==42)*k43 + ((combination%47)==41)*k42 + ((combination%47)==40)*k41 + ((combination%47)==39)*k40 +
1524 ((combination%47)==38)*k39 + ((combination%47)==37)*k38 + ((combination%47)==36)*k37 + ((combination%47)==35)*k36 + ((combination%47)==34)*k35 + ((combination%47)==33)*k34 + ((combination%47)==32)*k33 + ((combination%47)==31)*k32 + ((combination%47)==30)*k31 + ((combination%47)==29)*k30 +
1525 ((combination%47)==28)*k29 + ((combination%47)==27)*k28 + ((combination%47)==26)*k27 + ((combination%47)==25)*k26 + ((combination%47)==24)*k25 + ((combination%47)==23)*k24 + ((combination%47)==22)*k23 + ((combination%47)==21)*k22 + ((combination%47)==20)*k21 + ((combination%47)==19)*k20 +
1526 ((combination%47)==18)*k19 + ((combination%47)==17)*k18 + ((combination%47)==16)*k17 + ((combination%47)==15)*k16 + ((combination%47)==14)*k15 + ((combination%47)==13)*k14 + ((combination%47)==12)*k13 + ((combination%47)==11)*k12 + ((combination%47)==10)*k11 + ((combination%47)==9)*k10 +
1527 ((combination%47)==8)*k9 + ((combination%47)==7)*k8 + ((combination%47)==6)*k7 + ((combination%47)==5)*k6 + ((combination%47)==4)*k5 + ((combination%47)==3)*k4 + ((combination%47)==2)*k3 + ((combination%47)==1)*k2 + ((combination%47)==0)*k1);
1528 for(l=0; l<Kcombo; l++) {
1529 // D1 has a different number of rows 'l' than D2 and D3 so you need to use something like 'switch' to check combination%3 without going out of bound of J1
1530 switch(combination%47) {
1532 D[j][Kcounter+l] = D1[l];
1535 D[j][Kcounter+l] = D2[l];
1538 D[j][Kcounter+l] = D3[l];
1541 D[j][Kcounter+l] = D4[l];
1544 D[j][Kcounter+l] = D5[l];
1547 D[j][Kcounter+l] = D6[l];
1550 D[j][Kcounter+l] = D7[l];
1553 D[j][Kcounter+l] = D8[l];
1556 D[j][Kcounter+l] = D9[l];
1559 D[j][Kcounter+l] = D10[l];
1562 D[j][Kcounter+l] = D11[l];
1565 D[j][Kcounter+l] = D12[l];
1568 D[j][Kcounter+l] = D13[l];
1571 D[j][Kcounter+l] = D14[l];
1574 D[j][Kcounter+l] = D15[l];
1577 D[j][Kcounter+l] = D16[l];
1580 D[j][Kcounter+l] = D17[l];
1583 D[j][Kcounter+l] = D18[l];
1586 D[j][Kcounter+l] = D19[l];
1589 D[j][Kcounter+l] = D20[l];
1592 D[j][Kcounter+l] = D21[l];
1595 D[j][Kcounter+l] = D22[l];
1598 D[j][Kcounter+l] = D23[l];
1601 D[j][Kcounter+l] = D24[l];
1604 D[j][Kcounter+l] = D25[l];
1607 D[j][Kcounter+l] = D26[l];
1610 D[j][Kcounter+l] = D27[l];
1613 D[j][Kcounter+l] = D28[l];
1616 D[j][Kcounter+l] = D29[l];
1619 D[j][Kcounter+l] = D30[l];
1622 D[j][Kcounter+l] = D31[l];
1625 D[j][Kcounter+l] = D32[l];
1628 D[j][Kcounter+l] = D33[l];
1631 D[j][Kcounter+l] = D34[l];
1634 D[j][Kcounter+l] = D35[l];
1637 D[j][Kcounter+l] = D36[l];
1640 D[j][Kcounter+l] = D37[l];
1643 D[j][Kcounter+l] = D38[l];
1646 D[j][Kcounter+l] = D39[l];
1649 D[j][Kcounter+l] = D40[l];
1652 D[j][Kcounter+l] = D41[l];
1655 D[j][Kcounter+l] = D42[l];
1658 D[j][Kcounter+l] = D43[l];
1661 D[j][Kcounter+l] = D44[l];
1664 D[j][Kcounter+l] = D45[l];
1667 D[j][Kcounter+l] = D46[l];
1670 D[j][Kcounter+l] = D47[l];
1676 for(m=0; m<Kcombo; m++) {
1677 // J1 has a different number of rows 'l' than J2 and J3 so you need to use something like 'switch' to check combination%3 without going out of bound of J1
1678 switch(combination%47) {
1680 J[j][Kcounter+l][Kcounter+m] = J1[l][m];
1683 J[j][Kcounter+l][Kcounter+m] = J2[l][m];
1686 J[j][Kcounter+l][Kcounter+m] = J3[l][m];
1689 J[j][Kcounter+l][Kcounter+m] = J4[l][m];
1692 J[j][Kcounter+l][Kcounter+m] = J5[l][m];
1695 J[j][Kcounter+l][Kcounter+m] = J6[l][m];
1698 J[j][Kcounter+l][Kcounter+m] = J7[l][m];
1701 J[j][Kcounter+l][Kcounter+m] = J8[l][m];
1704 J[j][Kcounter+l][Kcounter+m] = J9[l][m];
1707 J[j][Kcounter+l][Kcounter+m] = J10[l][m];
1710 J[j][Kcounter+l][Kcounter+m] = J11[l][m];
1713 J[j][Kcounter+l][Kcounter+m] = J12[l][m];
1716 J[j][Kcounter+l][Kcounter+m] = J13[l][m];
1719 J[j][Kcounter+l][Kcounter+m] = J14[l][m];
1722 J[j][Kcounter+l][Kcounter+m] = J15[l][m];
1725 J[j][Kcounter+l][Kcounter+m] = J16[l][m];
1728 J[j][Kcounter+l][Kcounter+m] = J17[l][m];
1731 J[j][Kcounter+l][Kcounter+m] = J18[l][m];
1734 J[j][Kcounter+l][Kcounter+m] = J19[l][m];
1737 J[j][Kcounter+l][Kcounter+m] = J20[l][m];
1740 J[j][Kcounter+l][Kcounter+m] = J21[l][m];
1743 J[j][Kcounter+l][Kcounter+m] = J22[l][m];
1746 J[j][Kcounter+l][Kcounter+m] = J23[l][m];
1749 J[j][Kcounter+l][Kcounter+m] = J24[l][m];
1752 J[j][Kcounter+l][Kcounter+m] = J25[l][m];
1755 J[j][Kcounter+l][Kcounter+m] = J26[l][m];
1758 J[j][Kcounter+l][Kcounter+m] = J27[l][m];
1761 J[j][Kcounter+l][Kcounter+m] = J28[l][m];
1764 J[j][Kcounter+l][Kcounter+m] = J29[l][m];
1767 J[j][Kcounter+l][Kcounter+m] = J30[l][m];
1770 J[j][Kcounter+l][Kcounter+m] = J31[l][m];
1773 J[j][Kcounter+l][Kcounter+m] = J32[l][m];
1776 J[j][Kcounter+l][Kcounter+m] = J33[l][m];
1779 J[j][Kcounter+l][Kcounter+m] = J34[l][m];
1782 J[j][Kcounter+l][Kcounter+m] = J35[l][m];
1785 J[j][Kcounter+l][Kcounter+m] = J36[l][m];
1788 J[j][Kcounter+l][Kcounter+m] = J37[l][m];
1791 J[j][Kcounter+l][Kcounter+m] = J38[l][m];
1794 J[j][Kcounter+l][Kcounter+m] = J39[l][m];
1797 J[j][Kcounter+l][Kcounter+m] = J40[l][m];
1800 J[j][Kcounter+l][Kcounter+m] = J41[l][m];
1803 J[j][Kcounter+l][Kcounter+m] = J42[l][m];
1806 J[j][Kcounter+l][Kcounter+m] = J43[l][m];
1809 J[j][Kcounter+l][Kcounter+m] = J44[l][m];
1812 J[j][Kcounter+l][Kcounter+m] = J45[l][m];
1815 J[j][Kcounter+l][Kcounter+m] = J46[l][m];
1818 J[j][Kcounter+l][Kcounter+m] = J47[l][m];
1827 for(l=0; l<n1; l++) { // assuming n1=n2=n3
1828 h[j][k*n1+l] = (((combination%47)==46)*h47[l] + ((combination%47)==45)*h46[l] + ((combination%47)==44)*h45[l] + ((combination%47)==43)*h44[l] + ((combination%47)==42)*h43[l] + ((combination%47)==41)*h42[l] + ((combination%47)==40)*h41[l] + ((combination%47)==39)*h40[l] +
1829 ((combination%47)==38)*h39[l] + ((combination%47)==37)*h38[l] + ((combination%47)==36)*h37[l] + ((combination%47)==35)*h36[l] + ((combination%47)==34)*h35[l] + ((combination%47)==33)*h34[l] + ((combination%47)==32)*h33[l] + ((combination%47)==31)*h32[l] + ((combination%47)==30)*h31[l] + ((combination%47)==29)*h30[l] +
1830 ((combination%47)==28)*h29[l] + ((combination%47)==27)*h28[l] + ((combination%47)==26)*h27[l] + ((combination%47)==25)*h26[l] + ((combination%47)==24)*h25[l] + ((combination%47)==23)*h24[l] + ((combination%47)==22)*h23[l] + ((combination%47)==21)*h22[l] + ((combination%47)==20)*h21[l] + ((combination%47)==19)*h20[l] +
1831 ((combination%47)==18)*h19[l] + ((combination%47)==17)*h18[l] + ((combination%47)==16)*h17[l] + ((combination%47)==15)*h16[l] + ((combination%47)==14)*h15[l] + ((combination%47)==13)*h14[l] + ((combination%47)==12)*h13[l] + ((combination%47)==11)*h12[l] + ((combination%47)==10)*h11[l] + ((combination%47)==9)*h10[l] +
1832 ((combination%47)==8)*h9[l] + ((combination%47)==7)*h8[l] + ((combination%47)==6)*h7[l] + ((combination%47)==5)*h6[l] + ((combination%47)==4)*h5[l] + ((combination%47)==3)*h4[l] + ((combination%47)==2)*h3[l] + ((combination%47)==1)*h2[l] + ((combination%47)==0)*h1[l]);
1834 // only filling the K[j] first rows of G and GBar here corresponding to the basis for D and J
1835 for(l=0; l<Kcombo; l++) {
1836 for(m=0; m<n1; m++) { // assuming n1=n2=n3
1837 G[j][Kcounter+l][k*n1+m] = (((combination%47)==46)*G47[l][m] + ((combination%47)==45)*G46[l][m] + ((combination%47)==44)*G45[l][m] + ((combination%47)==43)*G44[l][m] + ((combination%47)==42)*G43[l][m] + ((combination%47)==41)*G42[l][m] + ((combination%47)==40)*G41[l][m] + ((combination%47)==39)*G40[l][m] +
1838 ((combination%47)==38)*G39[l][m] + ((combination%47)==37)*G38[l][m] + ((combination%47)==36)*G37[l][m] + ((combination%47)==35)*G36[l][m] + ((combination%47)==34)*G35[l][m] + ((combination%47)==33)*G34[l][m] + ((combination%47)==32)*G33[l][m] + ((combination%47)==31)*G32[l][m] + ((combination%47)==30)*G31[l][m] + ((combination%47)==29)*G30[l][m] +
1839 ((combination%47)==28)*G29[l][m] + ((combination%47)==27)*G28[l][m] + ((combination%47)==26)*G27[l][m] + ((combination%47)==25)*G26[l][m] + ((combination%47)==24)*G25[l][m] + ((combination%47)==23)*G24[l][m] + ((combination%47)==22)*G23[l][m] + ((combination%47)==21)*G22[l][m] + ((combination%47)==20)*G21[l][m] + ((combination%47)==19)*G20[l][m] +
1840 ((combination%47)==18)*G19[l][m] + ((combination%47)==17)*G18[l][m] + ((combination%47)==16)*G17[l][m] + ((combination%47)==15)*G16[l][m] + ((combination%47)==14)*G15[l][m] + ((combination%47)==13)*G14[l][m] + ((combination%47)==12)*G13[l][m] + ((combination%47)==11)*G12[l][m] + ((combination%47)==10)*G11[l][m] + ((combination%47)==9)*G10[l][m] +
1841 ((combination%47)==8)*G9[l][m] + ((combination%47)==7)*G8[l][m] + ((combination%47)==6)*G7[l][m] + ((combination%47)==5)*G6[l][m] + ((combination%47)==4)*G5[l][m] + ((combination%47)==3)*G4[l][m] + ((combination%47)==2)*G3[l][m] + ((combination%47)==1)*G2[l][m] + ((combination%47)==0)*G1[l][m]);
1842 GBar[j][Kcounter+l][k*n1+m] = (((combination%47)==46)*GBar47[l][m] + ((combination%47)==45)*GBar46[l][m] + ((combination%47)==44)*GBar45[l][m] + ((combination%47)==43)*GBar44[l][m] + ((combination%47)==42)*GBar43[l][m] + ((combination%47)==41)*GBar42[l][m] + ((combination%47)==40)*GBar41[l][m] + ((combination%47)==39)*GBar40[l][m] +
1843 ((combination%47)==38)*GBar39[l][m] + ((combination%47)==37)*GBar38[l][m] + ((combination%47)==36)*GBar37[l][m] + ((combination%47)==35)*GBar36[l][m] + ((combination%47)==34)*GBar35[l][m] + ((combination%47)==33)*GBar34[l][m] + ((combination%47)==32)*GBar33[l][m] + ((combination%47)==31)*GBar32[l][m] + ((combination%47)==30)*GBar31[l][m] + ((combination%47)==29)*GBar30[l][m] +
1844 ((combination%47)==28)*GBar29[l][m] + ((combination%47)==27)*GBar28[l][m] + ((combination%47)==26)*GBar27[l][m] + ((combination%47)==25)*GBar26[l][m] + ((combination%47)==24)*GBar25[l][m] + ((combination%47)==23)*GBar24[l][m] + ((combination%47)==22)*GBar23[l][m] + ((combination%47)==21)*GBar22[l][m] + ((combination%47)==20)*GBar21[l][m] + ((combination%47)==19)*GBar20[l][m] +
1845 ((combination%47)==18)*GBar19[l][m] + ((combination%47)==17)*GBar18[l][m] + ((combination%47)==16)*GBar17[l][m] + ((combination%47)==15)*GBar16[l][m] + ((combination%47)==14)*GBar15[l][m] + ((combination%47)==13)*GBar14[l][m] + ((combination%47)==12)*GBar13[l][m] + ((combination%47)==11)*GBar12[l][m] + ((combination%47)==10)*GBar11[l][m] + ((combination%47)==9)*GBar10[l][m] +
1846 ((combination%47)==8)*GBar9[l][m] + ((combination%47)==7)*GBar8[l][m] + ((combination%47)==6)*GBar7[l][m] + ((combination%47)==5)*GBar6[l][m] + ((combination%47)==4)*GBar5[l][m] + ((combination%47)==3)*GBar4[l][m] + ((combination%47)==2)*GBar3[l][m] + ((combination%47)==1)*GBar2[l][m] + ((combination%47)==0)*GBar1[l][m]);
1849 Kcounter = Kcounter + Kcombo;
1851 /* printf("intermediate G[%d]:\n", j); */
1852 /* printMatrix(G[j], N, N); */
1853 /* printf("intermediate GBar[%d]:\n", j); */
1854 /* printMatrix(GBar[j], N, N); */
1855 //memcpy(origG[j][k], G[j][k], N*sizeof(int)); memcpy(origGBar[j][k], GBar[j][k], N*sizeof(int));
1857 //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
1859 combination /= 47; // shift to the right by one (in base-47 arithmetic)
1863 // now need to fill the N-Kcounter remaining rows of G and GBar that are outside the spanning basis states of D and J
1865 for(k=0; k<(N/12); k++) {
1866 Kcombo = (((combination%47)==46)*k47 + ((combination%47)==45)*k46 + ((combination%47)==44)*k45 + ((combination%47)==43)*k44 + ((combination%47)==42)*k43 + ((combination%47)==41)*k42 + ((combination%47)==40)*k41 + ((combination%47)==39)*k40 +
1867 ((combination%47)==38)*k39 + ((combination%47)==37)*k38 + ((combination%47)==36)*k37 + ((combination%47)==35)*k36 + ((combination%47)==34)*k35 + ((combination%47)==33)*k34 + ((combination%47)==32)*k33 + ((combination%47)==31)*k32 + ((combination%47)==30)*k31 + ((combination%47)==29)*k30 +
1868 ((combination%47)==28)*k29 + ((combination%47)==27)*k28 + ((combination%47)==26)*k27 + ((combination%47)==25)*k26 + ((combination%47)==24)*k25 + ((combination%47)==23)*k24 + ((combination%47)==22)*k23 + ((combination%47)==21)*k22 + ((combination%47)==20)*k21 + ((combination%47)==19)*k20 +
1869 ((combination%47)==18)*k19 + ((combination%47)==17)*k18 + ((combination%47)==16)*k17 + ((combination%47)==15)*k16 + ((combination%47)==14)*k15 + ((combination%47)==13)*k14 + ((combination%47)==12)*k13 + ((combination%47)==11)*k12 + ((combination%47)==10)*k11 + ((combination%47)==9)*k10 +
1870 ((combination%47)==8)*k9 + ((combination%47)==7)*k8 + ((combination%47)==6)*k7 + ((combination%47)==5)*k6 + ((combination%47)==4)*k5 + ((combination%47)==3)*k4 + ((combination%47)==2)*k3 + ((combination%47)==1)*k2 + ((combination%47)==0)*k1);
1871 //printf("Kcounter=%d\n", Kcounter);
1872 // G and GBar rows that are outside the first 'k' spanning basis states
1873 for(l=Kcombo; l<n1; l++) { // assuming n1=n2=n3
1874 //printf("l=%d\n", l);
1875 for(m=0; m<n1; m++) { // assuming n1=n2=n3
1876 /* printf("m=%d\n", m); */
1877 /* printf("Kcounter+l=%d\n", Kcounter+l); */
1878 /* printf("k*n1+m=%d\n", k*n1+m); */
1879 G[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%47)==46)*G47[l][m] + ((combination%47)==45)*G46[l][m] + ((combination%47)==44)*G45[l][m] + ((combination%47)==43)*G44[l][m] + ((combination%47)==42)*G43[l][m] + ((combination%47)==41)*G42[l][m] + ((combination%47)==40)*G41[l][m] + ((combination%47)==39)*G40[l][m] +
1880 ((combination%47)==38)*G39[l][m] + ((combination%47)==37)*G38[l][m] + ((combination%47)==36)*G37[l][m] + ((combination%47)==35)*G36[l][m] + ((combination%47)==34)*G35[l][m] + ((combination%47)==33)*G34[l][m] + ((combination%47)==32)*G33[l][m] + ((combination%47)==31)*G32[l][m] + ((combination%47)==30)*G31[l][m] + ((combination%47)==29)*G30[l][m] +
1881 ((combination%47)==28)*G29[l][m] + ((combination%47)==27)*G28[l][m] + ((combination%47)==26)*G27[l][m] + ((combination%47)==25)*G26[l][m] + ((combination%47)==24)*G25[l][m] + ((combination%47)==23)*G24[l][m] + ((combination%47)==22)*G23[l][m] + ((combination%47)==21)*G22[l][m] + ((combination%47)==20)*G21[l][m] + ((combination%47)==19)*G20[l][m] +
1882 ((combination%47)==18)*G19[l][m] + ((combination%47)==17)*G18[l][m] + ((combination%47)==16)*G17[l][m] + ((combination%47)==15)*G16[l][m] + ((combination%47)==14)*G15[l][m] + ((combination%47)==13)*G14[l][m] + ((combination%47)==12)*G13[l][m] + ((combination%47)==11)*G12[l][m] + ((combination%47)==10)*G11[l][m] + ((combination%47)==9)*G10[l][m] +
1883 ((combination%47)==8)*G9[l][m] + ((combination%47)==7)*G8[l][m] + ((combination%47)==6)*G7[l][m] + ((combination%47)==5)*G6[l][m] + ((combination%47)==4)*G5[l][m] + ((combination%47)==3)*G4[l][m] + ((combination%47)==2)*G3[l][m] + ((combination%47)==1)*G2[l][m] + ((combination%47)==0)*G1[l][m]);
1884 GBar[j][Kcounter+l-Kcombo][k*n1+m] = (((combination%47)==46)*GBar47[l][m] + ((combination%47)==45)*GBar46[l][m] + ((combination%47)==44)*GBar45[l][m] + ((combination%47)==43)*GBar44[l][m] + ((combination%47)==42)*GBar43[l][m] + ((combination%47)==41)*GBar42[l][m] + ((combination%47)==40)*GBar41[l][m] + ((combination%47)==39)*GBar40[l][m] +
1885 ((combination%47)==38)*GBar39[l][m] + ((combination%47)==37)*GBar38[l][m] + ((combination%47)==36)*GBar37[l][m] + ((combination%47)==35)*GBar36[l][m] + ((combination%47)==34)*GBar35[l][m] + ((combination%47)==33)*GBar34[l][m] + ((combination%47)==32)*GBar33[l][m] + ((combination%47)==31)*GBar32[l][m] + ((combination%47)==30)*GBar31[l][m] + ((combination%47)==29)*GBar30[l][m] +
1886 ((combination%47)==28)*GBar29[l][m] + ((combination%47)==27)*GBar28[l][m] + ((combination%47)==26)*GBar27[l][m] + ((combination%47)==25)*GBar26[l][m] + ((combination%47)==24)*GBar25[l][m] + ((combination%47)==23)*GBar24[l][m] + ((combination%47)==22)*GBar23[l][m] + ((combination%47)==21)*GBar22[l][m] + ((combination%47)==20)*GBar21[l][m] + ((combination%47)==19)*GBar20[l][m] +
1887 ((combination%47)==18)*GBar19[l][m] + ((combination%47)==17)*GBar18[l][m] + ((combination%47)==16)*GBar17[l][m] + ((combination%47)==15)*GBar16[l][m] + ((combination%47)==14)*GBar15[l][m] + ((combination%47)==13)*GBar14[l][m] + ((combination%47)==12)*GBar13[l][m] + ((combination%47)==11)*GBar12[l][m] + ((combination%47)==10)*GBar11[l][m] + ((combination%47)==9)*GBar10[l][m] +
1888 ((combination%47)==8)*GBar9[l][m] + ((combination%47)==7)*GBar8[l][m] + ((combination%47)==6)*GBar7[l][m] + ((combination%47)==5)*GBar6[l][m] + ((combination%47)==4)*GBar5[l][m] + ((combination%47)==3)*GBar4[l][m] + ((combination%47)==2)*GBar3[l][m] + ((combination%47)==1)*GBar2[l][m] + ((combination%47)==0)*GBar1[l][m]);
1891 Kcounter = Kcounter + (n1-Kcombo);
1893 /* printf("intermediate G[%d]:\n", j); */
1894 /* printMatrix(G[j], N, N); */
1895 /* printf("intermediate GBar[%d]:\n", j); */
1896 /* printMatrix(GBar[j], N, N); */
1901 /*printf("G[%d]:\n", j);
1902 printMatrix(G[j], N, N);
1903 printf("GBar[%d]:\n", j);
1904 printMatrix(GBar[j], N, N);
1906 printf("h[%d]:\n", j);
1907 printVector(h[j], N);
1909 printf("J[%d]:\n", j);
1910 printMatrix(J[j], K[j], K[j]);
1912 printf("D[%d]:\n", j);
1913 printVector(D[j], K[j]);
1915 printf("Q[%d]=%d\n", j, Q[j]);*/
1920 while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
1922 if((Paulicounter+1) > N) {
1923 printf("Error: Number of Paulis is greater than N!\n");
1927 // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
1929 for(i=0; i<N; i++) {
1930 if(delta[Paulicounter][i]){
1931 omega[Paulicounter] += 3; // -I = I^3
1932 beta[Paulicounter][i] = delta[Paulicounter][i];
1933 gamma[Paulicounter][i] = delta[Paulicounter][i];
1937 /*printf("*******\n");
1938 printf("*******\n");
1939 printf("omega=%d\n", omega);
1941 printVector(gamma, N);
1943 printVector(beta, N);
1944 printf("*******\n");
1945 printf("*******\n");*/
1947 //for(j=0; j<pow(47,N/12); j++) { // the kets
1949 /*printf("========\n");
1950 printf("before:\n");
1951 printf("K=%d\n", K[j]);
1953 printVector(h[j], N);
1954 printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
1956 printMatrix(G[j], N, N);
1958 printMatrix(GBar[j], N, N);
1959 printf("Q=%d\n", Q[j]);
1961 printVector(D[j], K[j]);
1963 printMatrix(J[j], K[j], K[j]);*/
1964 //Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
1965 /*printf("\nafter:\n");
1966 printf("K=%d\n", K[j]);
1968 printVector(h[j], N);
1969 printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
1971 printMatrix(G[j], N, N);
1973 printMatrix(GBar[j], N, N);
1974 printf("Q=%d\n", Q[j]);
1976 printVector(D[j], K[j]);
1978 printMatrix(J[j], K[j], K[j]);*/
1985 double complex amplitude = 0.0 + 0.0*I;
1986 for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
1987 //printf("i=%d\n", i);
1989 randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
1991 origGamma = 1.0 + 0.0*I;
1993 for(k=0; k<Paulicounter; k++) {
1994 origGamma *= measurepauli(N, &origK, origh, origG, origGBar, &origQ, &origD, &origJ, omega[k], gamma[k], beta[k]);
1995 //printf("k=%d\n", k);
1997 /*printf("origK=%d\n", origK);
1999 printMatrix(origG, N, N);
2000 printf("origGBar:\n");
2001 printMatrix(origGBar, N, N);
2003 printVector(origh, N);*/
2005 double complex stabstateaverage = 0.0 + 0.0*I;
2007 for(j=0; j<pow(47,N/12); j++) {
2008 //printf("j=%d\n", j);
2009 double complex newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK, origh, origG, origGBar, origQ, origD, origJ);
2010 stabstateaverage = stabstateaverage + origGamma*Gamma[j]*newamplitude;
2012 amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
2014 deallocate_mem(&origG, N);
2015 deallocate_mem(&origGBar, N);
2017 deallocate_mem(&origJ, origK);
2021 printf("amplitude:\n");
2022 if(creal(amplitude+ZEROTHRESHOLD)>0)
2023 printf("%.10lf %c %.10lf I\n", cabs(creal(amplitude)), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
2025 printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
2029 for(i=0; i<PdN; i++)
2036 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
2039 int newomega, newalpha, newbeta, newgamma, newdelta;
2042 if(scanf("%d", &newomega) != EOF) {
2044 for(i=0; i<numqubits; i++) {
2045 if(scanf("%d %d %d %d", &newalpha, &newbeta, &newgamma, &newdelta) == EOF) {
2046 printf("Error: Too few input coeffs!\n");
2049 if(newalpha+newbeta+newgamma+newdelta > 1) {
2050 printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
2053 alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;