strongsim12_relerr.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <complex.h>
   6 #include <time.h>
   7 #include "matrix.h"
   8 #include "exponentialsum.h"
   9 #include "shrinkstar.h"
  10 #include "extend.h"
  11 #include "measurepauli.h"
  12 #include "innerproduct.h"
  13 #include "randomstabilizerstate.h"
  14
  15 #define ZEROTHRESHOLD (0.00000001)
  16
  17 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits);
  18
  19 // order of matrix elements is [row][column]!!!
  20
  21 int main(int argc, char *argv[])
  22 {
  23
  24   if(argc != 2) {
  25     printf("strongsim12_rellerr argument: \"number of stabilizer state samples\"\n");
  26     exit(0);
  27   }
  28
  29   int NUMSTABSTATESAMPLES = atoi(argv[1]);        // number of stabilizer state samples
  30
  31   int N;                  // number of qubits
  32   scanf("%d", &N);
  33
  34   if(N%12 != 0) {
  35     printf("'N' needs to be a multiple of 12 for a k=12 tensor factor decomposition!\n");
  36     return 1;
  37   }
  38
  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)
  40   scanf("%d", &T);
  41
  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
  44   int Paulicounter = 0;
  45
  46   int i, j, k, l, m;
  47
  48   FILE *fp;
  49   char buff[255];
  50
  51   srand((unsigned)time(NULL)); // seeding the random number generator for randomstabilizerstate()
  52
  53   fp = fopen("Pd.txt", "r");
  54
  55   if(fscanf(fp, "%s", buff) == EOF) {
  56     printf("Error: Pd.txt should start with the number N of P(d) of values tabulated.");
  57     return 1;
  58   }
  59
  60   double** Pd;
  61
  62   int PdN = atoi(buff);
  63
  64   Pd = calloc(PdN, sizeof(double*));
  65   for(i=0; i<PdN; i++)
  66     Pd[i] = calloc(PdN+1, sizeof(double));
  67
  68   double tmp;
  69
  70   for(i=1; i<PdN; i++) {
  71     tmp = 0.0;
  72     for(j=0; j<=i; j++) {
  73       if(fscanf(fp, "%s", buff) == EOF) {
  74         printf("Error: expected more values tabulated for P(d) for N=%d", PdN);
  75         return 1;
  76       }
  77       Pd[i][j] = atof(buff);
  78       //printf("%e ", Pd[i][j]);
  79       tmp += Pd[i][j];
  80     }
  81     //printf("\n");
  82     //printf("total=%f\n", tmp);
  83   }
  84
  85
  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);
  93
  94   // b60
  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} };
  96   // b66
  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} };
  98   // e6
  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}  };
 100   // o6
 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}  };
 102   // k6
 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} };
 104   // phiprime
 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}  };
 106   // phidoubleprime
 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}  };
 108
 109   // b60 * b60
 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));
 118   }
 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);
 132   Q1 = (QQ1+QQ1);
 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);
 137
 138   // b66 * b60
 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));
 147   }
 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);
 161   Q2 = (QQ1+QQ2);
 162   appendVector(DD2, DD1, D2, kk2, kk1);
 163   appendBlockMatrix(JJ2, JJ1, J2, kk2, kk2, kk1, kk1);*/
 164
 165   // e6 * b60
 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));
 174   }
 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);
 188   Q2 = (QQ1+QQ3);
 189   appendVector(DD3, DD1, D2, kk3, kk1);
 190   appendBlockMatrix(JJ3, JJ1, J2, kk3, kk3, kk1, kk1);
 191
 192   // o6 * b60
 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));
 201   }
 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);
 215   Q3 = (QQ1+QQ4);
 216   appendVector(DD4, DD1, D3, kk4, kk1);
 217   appendBlockMatrix(JJ4, JJ1, J3, kk4, kk4, kk1, kk1);
 218
 219   // k6 * b60
 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));
 228   }
 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);
 242   Q4 = (QQ1+QQ5);
 243   appendVector(DD5, DD1, D4, kk5, kk1);
 244   appendBlockMatrix(JJ5, JJ1, J4, kk5, kk5, kk1, kk1);
 245
 246   // phiprime * b60
 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));
 255   }
 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);
 269   Q5 = (QQ1+QQ6);
 270   appendVector(DD6, DD1, D5, kk6, kk1);
 271   appendBlockMatrix(JJ6, JJ1, J5, kk6, kk6, kk1, kk1);
 272
 273   // phidprime * b60
 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));
 282   }
 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);
 296   Q6 = (QQ1+QQ7);
 297   appendVector(DD7, DD1, D6, kk7, kk1);
 298   appendBlockMatrix(JJ7, JJ1, J6, kk7, kk7, kk1, kk1);
 299
 300   // b60 * b66
 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));
 309   }
 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);
 323   Q8 = (QQ2+QQ1);
 324   appendVector(DD1, DD2, D8, kk1, kk2);
 325   appendBlockMatrix(JJ1, JJ2, J8, kk1, kk1, kk2, kk2);*/
 326
 327   // b66 * b66
 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));
 336   }
 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);
 350   Q7 = (QQ2+QQ2);
 351   appendVector(DD2, DD2, D7, kk2, kk2);
 352   appendBlockMatrix(JJ2, JJ2, J7, kk2, kk2, kk2, kk2);
 353
 354   // e6 * b66
 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));
 363   }
 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);
 377   Q8 = (QQ2+QQ3);
 378   appendVector(DD3, DD2, D8, kk3, kk2);
 379   appendBlockMatrix(JJ3, JJ2, J8, kk3, kk3, kk2, kk2);
 380
 381   // o6 * b66
 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));
 390   }
 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);
 404   Q9 = (QQ2+QQ4);
 405   appendVector(DD4, DD2, D9, kk4, kk2);
 406   appendBlockMatrix(JJ4, JJ2, J9, kk4, kk4, kk2, kk2);
 407
 408   // k6 * b66
 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));
 417   }
 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);
 431   Q10 = (QQ2+QQ5);
 432   appendVector(DD5, DD2, D10, kk5, kk2);
 433   appendBlockMatrix(JJ5, JJ2, J10, kk5, kk5, kk2, kk2);
 434
 435   // phiprime * b66
 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));
 444   }
 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);
 458   Q11 = (QQ2+QQ6);
 459   appendVector(DD6, DD2, D11, kk6, kk2);
 460   appendBlockMatrix(JJ6, JJ2, J11, kk6, kk6, kk2, kk2);
 461
 462   // phidprime * b66
 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));
 471   }
 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);
 485   Q12 = (QQ2+QQ7);
 486   appendVector(DD7, DD2, D12, kk7, kk2);
 487   appendBlockMatrix(JJ7, JJ2, J12, kk7, kk7, kk2, kk2);
 488
 489   // b60 * e6
 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));
 498   }
 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);
 512   Q13 = (QQ3+QQ1);
 513   appendVector(DD1, DD3, D13, kk1, kk3);
 514   appendBlockMatrix(JJ1, JJ3, J13, kk1, kk1, kk3, kk3);
 515
 516   // b66 * e6
 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));
 525   }
 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);
 539   Q14 = (QQ3+QQ2);
 540   appendVector(DD2, DD3, D14, kk2, kk3);
 541   appendBlockMatrix(JJ2, JJ3, J14, kk2, kk2, kk3, kk3);
 542
 543   // e6 * e6
 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));
 552   }
 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);
 566   Q15 = (QQ3+QQ3);
 567   appendVector(DD3, DD3, D15, kk3, kk3);
 568   appendBlockMatrix(JJ3, JJ3, J15, kk3, kk3, kk3, kk3);
 569
 570   // o6 * e6
 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));
 579   }
 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);
 593   Q18 = (QQ3+QQ4);
 594   appendVector(DD4, DD3, D18, kk4, kk3);
 595   appendBlockMatrix(JJ4, JJ3, J18, kk4, kk4, kk3, kk3);*/
 596
 597   // k6 * e6
 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));
 606   }
 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);
 620   Q16 = (QQ3+QQ5);
 621   appendVector(DD5, DD3, D16, kk5, kk3);
 622   appendBlockMatrix(JJ5, JJ3, J16, kk5, kk5, kk3, kk3);
 623
 624   // phiprime * e6
 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));
 633   }
 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);
 647   Q17 = (QQ3+QQ6);
 648   appendVector(DD6, DD3, D17, kk6, kk3);
 649   appendBlockMatrix(JJ6, JJ3, J17, kk6, kk6, kk3, kk3);
 650
 651   // phidprime * e6
 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));
 660   }
 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);
 674   Q18 = (QQ3+QQ7);
 675   appendVector(DD7, DD3, D18, kk7, kk3);
 676   appendBlockMatrix(JJ7, JJ3, J18, kk7, kk7, kk3, kk3);
 677
 678   // b60 * o6
 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));
 687   }
 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);
 701   Q19 = (QQ4+QQ1);
 702   appendVector(DD1, DD4, D19, kk1, kk4);
 703   appendBlockMatrix(JJ1, JJ4, J19, kk1, kk1, kk4, kk4);
 704
 705   // b66 * o6
 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));
 714   }
 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);
 728   Q20 = (QQ4+QQ2);
 729   appendVector(DD2, DD4, D20, kk2, kk4);
 730   appendBlockMatrix(JJ2, JJ4, J20, kk2, kk2, kk4, kk4);
 731
 732   // e6 * o6
 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));
 741   }
 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);
 755   Q24 = (QQ4+QQ3);
 756   appendVector(DD3, DD4, D24, kk3, kk4);
 757   appendBlockMatrix(JJ3, JJ4, J24, kk3, kk3, kk4, kk4);*/
 758
 759   // o6 * o6
 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));
 768   }
 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);
 782   Q21 = (QQ4+QQ4);
 783   appendVector(DD4, DD4, D21, kk4, kk4);
 784   appendBlockMatrix(JJ4, JJ4, J21, kk4, kk4, kk4, kk4);
 785
 786   // k6 * o6
 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));
 795   }
 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);
 809   Q22 = (QQ4+QQ5);
 810   appendVector(DD5, DD4, D22, kk5, kk4);
 811   appendBlockMatrix(JJ5, JJ4, J22, kk5, kk5, kk4, kk4);
 812
 813   // phiprime * o6
 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));
 822   }
 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);
 836   Q23 = (QQ4+QQ6);
 837   appendVector(DD6, DD4, D23, kk6, kk4);
 838   appendBlockMatrix(JJ6, JJ4, J23, kk6, kk6, kk4, kk4);
 839
 840   // phidprime * o6
 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));
 849   }
 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);
 863   Q24 = (QQ4+QQ7);
 864   appendVector(DD7, DD4, D24, kk7, kk4);
 865   appendBlockMatrix(JJ7, JJ4, J24, kk7, kk7, kk4, kk4);
 866
 867   // b60 * k6
 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));
 876   }
 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);
 890   Q25 = (QQ5+QQ1);
 891   appendVector(DD1, DD5, D25, kk1, kk5);
 892   appendBlockMatrix(JJ1, JJ5, J25, kk1, kk1, kk5, kk5);
 893
 894   // b66 * k6
 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));
 903   }
 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);
 917   Q26 = (QQ5+QQ2);
 918   appendVector(DD2, DD5, D26, kk2, kk5);
 919   appendBlockMatrix(JJ2, JJ5, J26, kk2, kk2, kk5, kk5);
 920
 921   // e6 * k6
 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));
 930   }
 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);
 944   Q27 = (QQ5+QQ3);
 945   appendVector(DD3, DD5, D27, kk3, kk5);
 946   appendBlockMatrix(JJ3, JJ5, J27, kk3, kk3, kk5, kk5);
 947
 948   // o6 * k6
 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));
 957   }
 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);
 971   Q28 = (QQ5+QQ4);
 972   appendVector(DD4, DD5, D28, kk4, kk5);
 973   appendBlockMatrix(JJ4, JJ5, J28, kk4, kk4, kk5, kk5);
 974
 975   // k6 * k6
 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));
 984   }
 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);
 998   Q29 = (QQ5+QQ5);
 999   appendVector(DD5, DD5, D29, kk5, kk5);
1000   appendBlockMatrix(JJ5, JJ5, J29, kk5, kk5, kk5, kk5);
1001
1002   // phiprime * k6
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));
1011   }
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);
1025   Q30 = (QQ5+QQ6);
1026   appendVector(DD6, DD5, D30, kk6, kk5);
1027   appendBlockMatrix(JJ6, JJ5, J30, kk6, kk6, kk5, kk5);
1028
1029   // phidprime * k6
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));
1038   }
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);
1052   Q31 = (QQ5+QQ7);
1053   appendVector(DD7, DD5, D31, kk7, kk5);
1054   appendBlockMatrix(JJ7, JJ5, J31, kk7, kk7, kk5, kk5);
1055
1056   // b60 * phiprime
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));
1065   }
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);
1079   Q32 = (QQ6+QQ1);
1080   appendVector(DD1, DD6, D32, kk1, kk6);
1081   appendBlockMatrix(JJ1, JJ6, J32, kk1, kk1, kk6, kk6);
1082
1083   // b66 * phiprime
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));
1092   }
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);
1106   Q33 = (QQ6+QQ2);
1107   appendVector(DD2, DD6, D33, kk2, kk6);
1108   appendBlockMatrix(JJ2, JJ6, J33, kk2, kk2, kk6, kk6);
1109
1110   // e6 * phiprime
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));
1119   }
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);
1133   Q34 = (QQ6+QQ3);
1134   appendVector(DD3, DD6, D34, kk3, kk6);
1135   appendBlockMatrix(JJ3, JJ6, J34, kk3, kk3, kk6, kk6);
1136
1137   // o6 * phiprime
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));
1146   }
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);
1160   Q35 = (QQ6+QQ4);
1161   appendVector(DD4, DD6, D35, kk4, kk6);
1162   appendBlockMatrix(JJ4, JJ6, J35, kk4, kk4, kk6, kk6);
1163
1164   // k6 * phiprime
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));
1173   }
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);
1187   Q36 = (QQ6+QQ5);
1188   appendVector(DD5, DD6, D36, kk5, kk6);
1189   appendBlockMatrix(JJ5, JJ6, J36, kk5, kk5, kk6, kk6);
1190
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));
1200   }
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);
1214   Q37 = (QQ6+QQ6);
1215   appendVector(DD6, DD6, D37, kk6, kk6);
1216   appendBlockMatrix(JJ6, JJ6, J37, kk6, kk6, kk6, kk6);
1217
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));
1227   }
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);
1241   Q38 = (QQ6+QQ7);
1242   appendVector(DD7, DD6, D38, kk7, kk6);
1243   appendBlockMatrix(JJ7, JJ6, J38, kk7, kk7, kk6, kk6);
1244
1245   // b60 * phidprime
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));
1254   }
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);
1268   Q39 = (QQ7+QQ1);
1269   appendVector(DD1, DD7, D39, kk1, kk7);
1270   appendBlockMatrix(JJ1, JJ7, J39, kk1, kk1, kk7, kk7);
1271
1272   // b66 * phidprime
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));
1281   }
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);
1295   Q40 = (QQ7+QQ2);
1296   appendVector(DD2, DD7, D40, kk2, kk7);
1297   appendBlockMatrix(JJ2, JJ7, J40, kk2, kk2, kk7, kk7);
1298
1299   // e6 * phidprime
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));
1308   }
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);
1322   Q41 = (QQ7+QQ3);
1323   appendVector(DD3, DD7, D41, kk3, kk7);
1324   appendBlockMatrix(JJ3, JJ7, J41, kk3, kk3, kk7, kk7);
1325
1326   // o6 * phidprime
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));
1335   }
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);
1349   Q42 = (QQ7+QQ4);
1350   appendVector(DD4, DD7, D42, kk4, kk7);
1351   appendBlockMatrix(JJ4, JJ7, J42, kk4, kk4, kk7, kk7);
1352
1353   // k6 * phidprime
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));
1362   }
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);
1376   Q43 = (QQ7+QQ5);
1377   appendVector(DD5, DD7, D43, kk5, kk7);
1378   appendBlockMatrix(JJ5, JJ7, J43, kk5, kk5, kk7, kk7);
1379
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));
1389   }
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);
1403   Q44 = (QQ7+QQ6);
1404   appendVector(DD6, DD7, D44, kk6, kk7);
1405   appendBlockMatrix(JJ6, JJ7, J44, kk6, kk6, kk7, kk7);
1406
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));
1416   }
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);
1430   Q45 = (QQ7+QQ7);
1431   appendVector(DD7, DD7, D45, kk7, kk7);
1432   appendBlockMatrix(JJ7, JJ7, J45, kk7, kk7, kk7, kk7);
1433
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};
1439   int Q46 = 4;
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}};
1442
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};
1448   int Q47 = 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}};
1451
1452
1453
1454
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*));
1459
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));
1461
1462   K = calloc(pow(47,N/12), sizeof(int));
1463
1464   int origK, origQ, *origD;
1465   int **origJ;
1466   int **origG, **origGBar;
1467   int *origh;
1468   double complex origGamma;
1469
1470   int combination; // a particular combination from the linear combo of stabilizer states making up the tensor factors multiplied together
1471
1472
1473   for(j=0; j<pow(47,N/12); j++) { // there will be 47^(N/12) combinations when using k=12 tensor factors
1474
1475     combination = j;
1476
1477     K[j] = 0.0;
1478
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);
1485       combination /= 47;
1486     }
1487     combination = j;
1488
1489     Gamma[j] = 1.0;
1490
1491     G[j] = calloc(N, sizeof(int*)); GBar[j] = calloc(N, sizeof(int*));
1492     h[j] = calloc(N, sizeof(int));
1493
1494     if(K[j] > 0) {
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));
1498     }
1499
1500     for(k=0; k<N; k++) {
1501       G[j][k] = calloc(N, sizeof(int)); GBar[j][k] = calloc(N, sizeof(int));
1502     }
1503
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++) {
1507
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);
1513
1514
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);
1522
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) {
1531         case 0:
1532           D[j][Kcounter+l] = D1[l];
1533           break;
1534         case 1:
1535           D[j][Kcounter+l] = D2[l];
1536           break;
1537         case 2:
1538           D[j][Kcounter+l] = D3[l];
1539           break;
1540         case 3:
1541           D[j][Kcounter+l] = D4[l];
1542             break;
1543         case 4:
1544           D[j][Kcounter+l] = D5[l];
1545           break;
1546         case 5:
1547           D[j][Kcounter+l] = D6[l];
1548           break;
1549         case 6:
1550           D[j][Kcounter+l] = D7[l];
1551           break;
1552         case 7:
1553           D[j][Kcounter+l] = D8[l];
1554           break;
1555         case 8:
1556           D[j][Kcounter+l] = D9[l];
1557           break;
1558         case 9:
1559           D[j][Kcounter+l] = D10[l];
1560           break;
1561         case 10:
1562           D[j][Kcounter+l] = D11[l];
1563           break;
1564         case 11:
1565           D[j][Kcounter+l] = D12[l];
1566           break;
1567         case 12:
1568           D[j][Kcounter+l] = D13[l];
1569           break;
1570         case 13:
1571           D[j][Kcounter+l] = D14[l];
1572             break;
1573         case 14:
1574           D[j][Kcounter+l] = D15[l];
1575           break;
1576         case 15:
1577           D[j][Kcounter+l] = D16[l];
1578           break;
1579         case 16:
1580           D[j][Kcounter+l] = D17[l];
1581           break;
1582         case 17:
1583           D[j][Kcounter+l] = D18[l];
1584           break;
1585         case 18:
1586           D[j][Kcounter+l] = D19[l];
1587           break;
1588         case 19:
1589           D[j][Kcounter+l] = D20[l];
1590           break;
1591         case 20:
1592           D[j][Kcounter+l] = D21[l];
1593           break;
1594         case 21:
1595           D[j][Kcounter+l] = D22[l];
1596           break;
1597         case 22:
1598           D[j][Kcounter+l] = D23[l];
1599           break;
1600         case 23:
1601           D[j][Kcounter+l] = D24[l];
1602             break;
1603         case 24:
1604           D[j][Kcounter+l] = D25[l];
1605           break;
1606         case 25:
1607           D[j][Kcounter+l] = D26[l];
1608           break;
1609         case 26:
1610           D[j][Kcounter+l] = D27[l];
1611           break;
1612         case 27:
1613           D[j][Kcounter+l] = D28[l];
1614           break;
1615         case 28:
1616           D[j][Kcounter+l] = D29[l];
1617           break;
1618         case 29:
1619           D[j][Kcounter+l] = D30[l];
1620           break;
1621         case 30:
1622           D[j][Kcounter+l] = D31[l];
1623           break;
1624         case 31:
1625           D[j][Kcounter+l] = D32[l];
1626           break;
1627         case 32:
1628           D[j][Kcounter+l] = D33[l];
1629           break;
1630         case 33:
1631           D[j][Kcounter+l] = D34[l];
1632             break;
1633         case 34:
1634           D[j][Kcounter+l] = D35[l];
1635           break;
1636         case 35:
1637           D[j][Kcounter+l] = D36[l];
1638           break;
1639         case 36:
1640           D[j][Kcounter+l] = D37[l];
1641           break;
1642         case 37:
1643           D[j][Kcounter+l] = D38[l];
1644           break;
1645         case 38:
1646           D[j][Kcounter+l] = D39[l];
1647           break;
1648         case 39:
1649           D[j][Kcounter+l] = D40[l];
1650           break;
1651         case 40:
1652           D[j][Kcounter+l] = D41[l];
1653           break;
1654         case 41:
1655           D[j][Kcounter+l] = D42[l];
1656           break;
1657         case 42:
1658           D[j][Kcounter+l] = D43[l];
1659           break;
1660         case 43:
1661           D[j][Kcounter+l] = D44[l];
1662             break;
1663         case 44:
1664           D[j][Kcounter+l] = D45[l];
1665           break;
1666         case 45:
1667           D[j][Kcounter+l] = D46[l];
1668           break;
1669         case 46:
1670           D[j][Kcounter+l] = D47[l];
1671           break;
1672         default:
1673           printf("error");
1674           return 1;
1675           }
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) {
1679           case 0:
1680             J[j][Kcounter+l][Kcounter+m] = J1[l][m];
1681             break;
1682           case 1:
1683             J[j][Kcounter+l][Kcounter+m] = J2[l][m];
1684             break;
1685           case 2:
1686             J[j][Kcounter+l][Kcounter+m] = J3[l][m];
1687             break;
1688           case 3:
1689             J[j][Kcounter+l][Kcounter+m] = J4[l][m];
1690             break;
1691           case 4:
1692             J[j][Kcounter+l][Kcounter+m] = J5[l][m];
1693             break;
1694           case 5:
1695             J[j][Kcounter+l][Kcounter+m] = J6[l][m];
1696             break;
1697           case 6:
1698             J[j][Kcounter+l][Kcounter+m] = J7[l][m];
1699             break;
1700           case 7:
1701             J[j][Kcounter+l][Kcounter+m] = J8[l][m];
1702             break;
1703           case 8:
1704             J[j][Kcounter+l][Kcounter+m] = J9[l][m];
1705             break;
1706           case 9:
1707             J[j][Kcounter+l][Kcounter+m] = J10[l][m];
1708             break;
1709           case 10:
1710             J[j][Kcounter+l][Kcounter+m] = J11[l][m];
1711             break;
1712           case 11:
1713             J[j][Kcounter+l][Kcounter+m] = J12[l][m];
1714             break;
1715           case 12:
1716             J[j][Kcounter+l][Kcounter+m] = J13[l][m];
1717             break;
1718           case 13:
1719             J[j][Kcounter+l][Kcounter+m] = J14[l][m];
1720             break;
1721           case 14:
1722             J[j][Kcounter+l][Kcounter+m] = J15[l][m];
1723             break;
1724           case 15:
1725             J[j][Kcounter+l][Kcounter+m] = J16[l][m];
1726             break;
1727           case 16:
1728             J[j][Kcounter+l][Kcounter+m] = J17[l][m];
1729             break;
1730           case 17:
1731             J[j][Kcounter+l][Kcounter+m] = J18[l][m];
1732             break;
1733           case 18:
1734             J[j][Kcounter+l][Kcounter+m] = J19[l][m];
1735             break;
1736           case 19:
1737             J[j][Kcounter+l][Kcounter+m] = J20[l][m];
1738             break;
1739           case 20:
1740             J[j][Kcounter+l][Kcounter+m] = J21[l][m];
1741             break;
1742           case 21:
1743             J[j][Kcounter+l][Kcounter+m] = J22[l][m];
1744             break;
1745           case 22:
1746             J[j][Kcounter+l][Kcounter+m] = J23[l][m];
1747             break;
1748           case 23:
1749             J[j][Kcounter+l][Kcounter+m] = J24[l][m];
1750             break;
1751           case 24:
1752             J[j][Kcounter+l][Kcounter+m] = J25[l][m];
1753             break;
1754           case 25:
1755             J[j][Kcounter+l][Kcounter+m] = J26[l][m];
1756             break;
1757           case 26:
1758             J[j][Kcounter+l][Kcounter+m] = J27[l][m];
1759             break;
1760           case 27:
1761             J[j][Kcounter+l][Kcounter+m] = J28[l][m];
1762             break;
1763           case 28:
1764             J[j][Kcounter+l][Kcounter+m] = J29[l][m];
1765             break;
1766           case 29:
1767             J[j][Kcounter+l][Kcounter+m] = J30[l][m];
1768             break;
1769           case 30:
1770             J[j][Kcounter+l][Kcounter+m] = J31[l][m];
1771             break;
1772           case 31:
1773             J[j][Kcounter+l][Kcounter+m] = J32[l][m];
1774             break;
1775           case 32:
1776             J[j][Kcounter+l][Kcounter+m] = J33[l][m];
1777             break;
1778           case 33:
1779             J[j][Kcounter+l][Kcounter+m] = J34[l][m];
1780             break;
1781           case 34:
1782             J[j][Kcounter+l][Kcounter+m] = J35[l][m];
1783             break;
1784           case 35:
1785             J[j][Kcounter+l][Kcounter+m] = J36[l][m];
1786             break;
1787           case 36:
1788             J[j][Kcounter+l][Kcounter+m] = J37[l][m];
1789             break;
1790           case 37:
1791             J[j][Kcounter+l][Kcounter+m] = J38[l][m];
1792             break;
1793           case 38:
1794             J[j][Kcounter+l][Kcounter+m] = J39[l][m];
1795             break;
1796           case 39:
1797             J[j][Kcounter+l][Kcounter+m] = J40[l][m];
1798             break;
1799           case 40:
1800             J[j][Kcounter+l][Kcounter+m] = J41[l][m];
1801             break;
1802           case 41:
1803             J[j][Kcounter+l][Kcounter+m] = J42[l][m];
1804             break;
1805           case 42:
1806             J[j][Kcounter+l][Kcounter+m] = J43[l][m];
1807             break;
1808           case 43:
1809             J[j][Kcounter+l][Kcounter+m] = J44[l][m];
1810             break;
1811           case 44:
1812             J[j][Kcounter+l][Kcounter+m] = J45[l][m];
1813             break;
1814           case 45:
1815             J[j][Kcounter+l][Kcounter+m] = J46[l][m];
1816             break;
1817           case 46:
1818             J[j][Kcounter+l][Kcounter+m] = J47[l][m];
1819             break;
1820           default:
1821             printf("error");
1822             return 1;
1823           }
1824         }
1825       }
1826
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]);
1833       }
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]);
1847         }
1848       }
1849       Kcounter = Kcounter + Kcombo;
1850
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));
1856
1857       //memcpy(origJ[j][k], J[j][k], K[j]*sizeof(int));
1858
1859       combination /= 47; // shift to the right by one (in base-47 arithmetic)
1860     }
1861     //printf("!\n");
1862
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
1864     combination = 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]);
1889         }
1890       }
1891       Kcounter = Kcounter + (n1-Kcombo);
1892
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); */
1897
1898       combination /= 47;
1899     }
1900
1901     /*printf("G[%d]:\n", j);
1902     printMatrix(G[j], N, N);
1903     printf("GBar[%d]:\n", j);
1904     printMatrix(GBar[j], N, N);
1905
1906     printf("h[%d]:\n", j);
1907     printVector(h[j], N);
1908
1909     printf("J[%d]:\n", j);
1910     printMatrix(J[j], K[j], K[j]);
1911
1912     printf("D[%d]:\n", j);
1913     printVector(D[j], K[j]);
1914
1915     printf("Q[%d]=%d\n", j, Q[j]);*/
1916
1917   }
1918   //exit(0);
1919
1920   while(readPaulicoeffs(&omega[Paulicounter], alpha[Paulicounter], beta[Paulicounter], gamma[Paulicounter], delta[Paulicounter], N)) {
1921
1922     if((Paulicounter+1) > N) {
1923       printf("Error: Number of Paulis is greater than N!\n");
1924       return 1;
1925     }
1926
1927     // Let's break up the Ys into Xs and Zs in the order Z X, as required to pass to measurepauli()
1928     // Y_i = -I*Z*X
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];
1934       }
1935     }
1936
1937     /*printf("*******\n");
1938     printf("*******\n");
1939     printf("omega=%d\n", omega);
1940     printf("X:\n");
1941     printVector(gamma, N);
1942     printf("Z:\n");
1943     printVector(beta, N);
1944     printf("*******\n");
1945     printf("*******\n");*/
1946
1947     //for(j=0; j<pow(47,N/12); j++) { // the kets
1948
1949       /*printf("========\n");
1950       printf("before:\n");
1951       printf("K=%d\n", K[j]);
1952       printf("h:\n");
1953       printVector(h[j], N);
1954       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
1955       printf("G:\n");
1956       printMatrix(G[j], N, N);
1957       printf("GBar:\n");
1958       printMatrix(GBar[j], N, N);
1959       printf("Q=%d\n", Q[j]);
1960       printf("D:\n");
1961       printVector(D[j], K[j]);
1962       printf("J:\n");
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]);
1967       printf("h:\n");
1968       printVector(h[j], N);
1969       printf("Gamma[%d]=%lf+%lf\n", j, creal(Gamma[j]), cimag(Gamma[j]));
1970       printf("G:\n");
1971       printMatrix(G[j], N, N);
1972       printf("GBar:\n");
1973       printMatrix(GBar[j], N, N);
1974       printf("Q=%d\n", Q[j]);
1975       printf("D:\n");
1976       printVector(D[j], K[j]);
1977       printf("J:\n");
1978       printMatrix(J[j], K[j], K[j]);*/
1979
1980     //}
1981
1982     Paulicounter++;
1983   }
1984
1985   double complex amplitude = 0.0 + 0.0*I;
1986   for(i=0; i<NUMSTABSTATESAMPLES; i++) { // the bras
1987     //printf("i=%d\n", i);
1988
1989     randomstabilizerstate(N, &origK, &origh, &origG, &origGBar, &origQ, &origD, &origJ, Pd);
1990
1991     origGamma = 1.0 + 0.0*I;
1992
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);
1996   }
1997     /*printf("origK=%d\n", origK);
1998     printf("origG:\n");
1999     printMatrix(origG, N, N);
2000     printf("origGBar:\n");
2001     printMatrix(origGBar, N, N);
2002     printf("origh:\n");
2003     printVector(origh, N);*/
2004
2005     double complex stabstateaverage = 0.0 + 0.0*I;
2006
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;
2011     }
2012     amplitude = amplitude + conj(stabstateaverage)*stabstateaverage/((double)(NUMSTABSTATESAMPLES))*pow(2.0,T);
2013
2014     deallocate_mem(&origG, N);
2015     deallocate_mem(&origGBar, N);
2016     free(origh);
2017     deallocate_mem(&origJ, origK);
2018     free(origD);
2019   }
2020
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)));
2024   else
2025     printf("%.10lf %c %.10lf I\n", creal(amplitude), cimag(amplitude+ZEROTHRESHOLD)>0?'+':'-' , cabs(cimag(amplitude)));
2026
2027
2028
2029   for(i=0; i<PdN; i++)
2030     free(Pd[i]);
2031   free(Pd);
2032
2033   return 0;
2034 }
2035
2036 int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
2037 {
2038
2039   int newomega, newalpha, newbeta, newgamma, newdelta;
2040   int i;
2041
2042   if(scanf("%d", &newomega) != EOF) {
2043     *omega = newomega;
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");
2047         exit(0);
2048       }
2049       if(newalpha+newbeta+newgamma+newdelta > 1) {
2050         printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
2051         exit(0);
2052       }
2053       alpha[i] = newalpha; beta[i] = newbeta; gamma[i] = newgamma; delta[i] = newdelta;
2054     }
2055     return 1;
2056   } else
2057     return 0;
2058
2059 }