product_lib.cpp

Go to the documentation of this file.

//outdated, not in use
//Weilei Zeng Nov 21, 2018
//to implement quantum concatenated codes. There are several ways of concatenation, see hypergraph_product_code.pdf
 
 
#include "dist.h"
#include <itpp/itbase.h>
#include <itpp/itcomm.h>
#include <stdio.h>
//#include "my_lib.h"
//#include "concatenation_lib.h"
#include "product_lib.h"
//using namespace common;
 
 
  //constructor
 
ClassicalCode::ClassicalCode(){
}
ClassicalCode::ClassicalCode(itpp::GF2mat G, itpp::GF2mat H){
  G=G;H=H;
  return;
}
 
//distance estimator
int ClassicalCode::dist(){
  return rand_dist();
}
int ClassicalCode::min_weight_dist(){
  return common::min_wt_decoding(G);
}
int ClassicalCode::rand_dist(){
  return common::rand_dist(G);
}
 
 
 
//function
void ClassicalCode::info(){
  std::cout<<"Classical code: n = "<<n<<std::endl
       <<"codeword generating matrix G"<<G
       <<"parity check matrix H"<<H<<std::endl;
  return;
}
ClassicalCode ClassicalCode::dual(){
  ClassicalCode dual_code(H,G);
  return dual_code;
}
 
void ClassicalCode::full_rank(){
  G = common::make_it_full_rank(G);
  H = common::make_it_full_rank(H);
}
 
//generate sample code
void ClassicalCode::get_repetition_code(){
  H =  common::get_check_rept(n);
  G = itpp::GF2mat(itpp::ones_b(n), false);
  return;
}
void ClassicalCode::get_743_code(){
  H =  common::get_check_code743(n);
  G =  common::get_check_code734(n);
  return;
}
void ClassicalCode::get_734_code(){
  H =  common::get_check_code734(n);
  G =  common::get_check_code743(n);
  return;
}
 
 
 
CSSCode::CSSCode(){}
CSSCode::CSSCode(int na, int Gax_row, int id_Gax, int Gaz_row, int id_Gaz){
    n=na;Gx_row=Gax_row; id_Gx=id_Gax;
    Gz_row=Gaz_row;id_Gz=id_Gaz;
    is_defined=1;
}
int CSSCode::generate_by_id(int debug){
    return generate_code(Gx, Gz, n, Gx_row, id_Gx, Gz_row, id_Gz, debug);
}
 
int CSSCode::getRandomCode(){
  return getRandomQuantumCode(n, Gx_row, Gz_row, Gx, Gz, Cx, Cz);
}
 
int CSSCode::getGoodCode(int debug){
  return getGoodQuantumCode(n, Gx_row, Gz_row, Gx, Gz, Cx, Cz, debug);
}
 
bool CSSCode::is_valid(){
  if ( is_C_defined ) {
    return common::is_quantum_code(Gx, Gz, Cx, Cz);
  }
  return common::is_quantum_code(Gx, Gz);
}
 
void CSSCode::dist(){
  dx = rand_dist_x();
  dz = rand_dist_z();
  return;
}
 
int CSSCode::min_weight_dist_x(){
  return  common::min_wt_decoding(Cx, Gx);
}
int CSSCode::min_weight_dist_z(){
  return  common::min_wt_decoding(Cz, Gz);
}
 
int CSSCode::rand_dist_x(){
  return common::quantum_dist_v2(Gx, Gz);
}
int CSSCode::rand_dist_z(){
  return common::quantum_dist_v2(Gx, Gz, 1);
}
 
void CSSCode::get_713_code(){
  Gx=common::get_check_code743(n);
  Gz=common::get_check_code743(n);
  return;
 
}
 
 
itpp::GF2mat remove_col(itpp::GF2mat G, int col){
  int n = G.cols();
  if ( col ==0 ) return G.get_submatrix(0,1,G.rows()-1,G.cols()-1);
  if ( col == n-1 ) return G.get_submatrix(0,0,G.rows()-1,G.cols()-2);
  return G.get_submatrix(0,0,G.rows()-1,col-1).concatenate_horizontal(
                                      G.get_submatrix(0,col+1,G.rows()-1,G.cols()-1)
                                      );  
}
 
void remove_singleton(itpp::GF2mat &Gx, itpp::GF2mat &Gz){
  //remove zero columns in Gx and Gz
  //not in use, just discard code with distance 1, easier solution
  int n = Gx.cols();
  itpp::bvec to_remove(n);//1 fro remove, 0 remain  
  for ( int i=0;i<n;i++){
    if (itpp::GF2mat(Gx.get_col(i)).is_zero()) to_remove.set(i,1);
    if (itpp::GF2mat(Gz.get_col(i)).is_zero()) to_remove.set(i,1);
  }
  for ( int i=0;i<n;i++){
    if ( to_remove(n-i-1) ){
    Gx=remove_col(Gx,n-i-1);
    Gz=remove_col(Gz,n-i-1);
      }
  }
  Gx=common::make_it_full_rank(Gx);
  Gz=common::make_it_full_rank(Gz);
  std::cout<<"singleton removed"<<std::endl;
  return;
}
 
int getRandomQuantumCode(int n,int Gx_row,int Gz_row, itpp::GF2mat &Gx,itpp::GF2mat &Gz, itpp::GF2mat &Cx,itpp::GF2mat &Cz){
 
  Gx = itpp::GF2mat(Gx_row,n);
  Gz = itpp::GF2mat(Gz_row,n);
  for ( int i =0;i<Gx_row;i++){//random G_x
    Gx.set_row(i,itpp::randb(n));//equally 0 and 1s
  }
  //Gx might not be full rank at this point
  
  itpp::GF2mat T,U; itpp::ivec P;
  int rank_of_Gx = Gx.transpose().T_fact(T,U,P);
  itpp::GF2mat Q=T.get_submatrix(rank_of_Gx,0,n-1,n-1);
 
 
  itpp::GF2mat alpha(Gz_row,Q.rows()); //a random binary matrix to select G_z
  for ( int i=0;i<Gz_row;i++){
    alpha.set_row(i,itpp::randb(Q.rows()));
  }
  Gz=alpha*Q;
  Cx=common::getC(Gx,Gz);
  Cz=common::getC(Gx,Gz,1);
  //  if (! is_quantum_code(Gx,Gz,Cx,Cz)) throw "invalid code";
 
  return 0;
}
 
int getGoodQuantumCode(int n,int Gx_row,int Gz_row, itpp::GF2mat &Gx,itpp::GF2mat &Gz, itpp::GF2mat &Cx,itpp::GF2mat &Cz, int debug){
  // return best codes among multip trial
  //repeat multiple times to get the best distance
  itpp::GF2mat Gx_temp, Gz_temp,Cx_temp,Cz_temp;
  int search_trial=1000;
  int flag_find_good_code=0;
  for ( int i =0; i<search_trial; i++){
 
    getRandomQuantumCode( n, Gx_row,Gz_row, Gx_temp, Gz_temp,Cx_temp,Cz_temp);
    //check distance and update if get larger distance
    int dx = common::quantum_dist_v2(Gx_temp,Gz_temp);
    if ( dx >1 ){
      int dz = common::quantum_dist_v2(Gx_temp,Gz_temp,1);
      if (dz >1 ){
    flag_find_good_code=1;
    //  Gx = Gx_temp; Gz = Gz_temp; Cx = Cx_temp; Cz = Cz_temp;
    if (debug) std::cout<<common::blue_text("get good code when i =")<<i<<std::endl;
    break;
    //  return 0;
      }
    }
  }
  Gx = Gx_temp; Gz = Gz_temp; Cx = Cx_temp; Cz = Cz_temp;
 
  //  if ( flag_find_good_code){
  if ( debug ) std::cout<<"Gx 1st row:"<<Gx.get_row(0)<<std::endl; // for debug the random seed
  //make sure both Gx and Gz are full rank
  if ( Gx.row_rank() < Gx.rows() ) {
    if (debug ) std::cout<<"getGoodQuantumCode: Gx not full rank. now make it full rank"<<std::endl;
    Gx = common::make_it_full_rank(Gx);
  }  
  if ( Gz.row_rank() < Gz.rows() ) {
    if (debug) std::cout<<"getGoodQuantumCode: Gz not full rank. now make it full rank"<<std::endl;
    Gz = common::make_it_full_rank(Gz);
  }
 
  if ( debug) if ( ! flag_find_good_code ) std::cout<<common::color_text("didn't find good code after ")<<search_trial<<" trials"<<std::endl;
 
 
  return 0;
}
 
void set_submatrix(itpp::GF2mat & G, itpp::GF2mat sub, int row, int col){
  //put sub into G, start from (row,col)
  for ( int i =0 ; i < sub.rows(); i ++)
    for ( int j = 0; j< sub.cols(); j++ ){
      G.set(i+row, j+col, sub.get(i,j));
    }
  return;
}
 
 
int is_row_reduced_echelon_form(itpp::GF2mat & alpha_Gaz, int debug = 0){
  //check it column by column, from bottom to top
  int get_one=0; //flag on if hit one in that column
  int position_one=-1;//position for one in that column
  int columns_one=0;//columns has a single one. exit for loop when reach alpha_Gaz.rows()
  for ( int i = 0; i<alpha_Gaz.cols();i++){
    get_one=0;
    for ( int j = alpha_Gaz.rows()-1; j > -1; j--){
      if (alpha_Gaz.get(j,i)){
    if (get_one){
      //get one twice in that column
      if (debug) std::cout<<"get one twice in that column i="<<i<<std::endl;
      //      std::cout<<"*";
      return 0;
    }else{
      if (j <= position_one){
        //skip this column, this column is not independent
        if (debug) std::cout<<"break the inner for loop for dependent column i = "<<i<<std::endl;
        break;
        //      continue;
      }else {
        get_one=1;
        position_one=j;
        columns_one++;
      }
    }
      }
    }
    //if (debug) std::cout<<"broke the inner for loop"<<std::endl;
    if ( columns_one == alpha_Gaz.rows() ){
      break;
    }
  }
  if ( columns_one < alpha_Gaz.rows() ){
    if (debug) std::cout<<"columns_one:"<<columns_one<<" is not full rank"<<alpha_Gaz.rows()<<std::endl;
    //    std::cout<<"*";
    return 0;
  }
  return 1;
}
 
// generate all code with size na systematically
int generate_code(itpp::GF2mat & Gax, itpp::GF2mat & Gaz, int na, int Gax_row, int id_Gax, int Gaz_row, int id_Gaz, int debug){
  if (debug) std::cout<<na<<","<<Gax_row<<","<<Gaz_row<<std::endl;
  //sanity check
  if (Gaz_row+Gax_row > na-1){
    std::cout<<"generate_code: no logical qubit"<<std::endl;
    throw 2;
  }
  const int id_Gax_MAX = (int) pow(2,  Gax_row * (na-Gax_row) ) -1 ; //maximun all one
  if ( id_Gax <1 || id_Gax > id_Gax_MAX ) {
    std::cout<<"illegal id_Gax: "<<id_Gax<<", id_Gax_MAX = "<<id_Gax_MAX<<std::endl;
    throw 2;
  }
  const int id_Gaz_MAX = (int) pow(2, Gaz_row*(na - Gax_row)) - 1; //maximun all one
  if ( id_Gaz < 1 || id_Gaz > id_Gaz_MAX ){
    std::cout<<"illegal id_Gaz: "<<id_Gaz<<", id_Gaz_MAX = "<<id_Gaz_MAX<<std::endl;
    throw 2;
  }
  //remove duplicate cases for id_Gax and id_Gaz
 
  itpp::GF2mat beta_Gaz = itpp::GF2mat(itpp::dec2bin(Gaz_row*(na-Gax_row),id_Gaz),false);
  itpp::GF2mat alpha_Gaz(Gaz_row, na-Gax_row);
  if ( debug ) std::cout<<"beta_Gaz = "<<beta_Gaz<<std::endl;
  for ( int i =0;i<Gaz_row;i++){
    if (debug) std::cout<<"set submatrix i = "<<i<<std::endl<<beta_Gaz.get_submatrix(0,i*(na-Gax_row),0, (i+1)*(na-Gax_row)-1)<<std::endl;
    set_submatrix(alpha_Gaz, beta_Gaz.get_submatrix(0,i*(na-Gax_row),0, (i+1)*(na-Gax_row)-1), i,0);
  } 
  if (debug) std::cout<<"alpha_Gaz"<<alpha_Gaz<<std::endl;
  /* decreasing order is ensured in reduced row echelon form, hence not checked here
  for ( int i =0;i<Gaz_row-1;i++){
    if ( itpp::bin2dec(alpha_Gaz.get_row(i)) <= itpp::bin2dec(alpha_Gaz.get_row(i+1))){
      if (debug) std::cout<< "duplicate Gaz with this id_Gaz. no calculation needed. alpha_Gaz/id_Gaz must be in decreasing order"<<std::endl;
      return 2;
    }
    }*/
  //make sure alpa_Gaz is in reduce row echelon form, to remove duplicate cases. return 2 if not in the form
  //this duplicate the check to make sure alpha_Gaz is in decreasing order
  if ( ! is_row_reduced_echelon_form( alpha_Gaz, debug) ) return 2;
 
 
  //finish check
 
 
 
  Gax = itpp::GF2mat(Gax_row,na);
  // identity matrix in the left part to make it reduce row echelon form.
  set_submatrix(Gax,itpp::gf2dense_eye(Gax_row),0,0);
  //  if (debug) std::cout<<"Gax"<<Gax<<std::endl;
 
 
  itpp::GF2mat alpha_Gax = itpp::GF2mat( itpp::dec2bin(Gax_row*(na-Gax_row), id_Gax), false);//false for row vector
  if (debug) std::cout<<"alpha_Gax give the right part of Gax"<<std::endl<<alpha_Gax<<std::endl;
  for ( int i = 0 ; i < Gax_row; i++){
    set_submatrix(Gax,alpha_Gax.get_submatrix(0, i*(na-Gax_row), 0, (i+1)*(na-Gax_row)-1), i, Gax_row);
  }
  if (debug) std::cout<<"Gax"<<Gax<<std::endl;
 
 
  //remove duplicate in id_Gax. They could be equal, but permute any two rows give equivalent code, so enfore all rows ( in the right part ) in decreasing order
  for ( int i =0;i<Gax_row-1;i++){
    if ( itpp::bin2dec(Gax.get_submatrix(0,Gax_row,Gax_row-1,na-1).get_row(i)) 
     < itpp::bin2dec(Gax.get_submatrix(0,Gax_row,Gax_row-1,na-1).get_row(i+1)) ){
      //itpp::bin2dec(alpha_Gaz.get_row(i+1))){
      if (debug) std::cout<< "duplicate Gax with this id_Gax. no calculation needed. id_Gax must be in decreasing order. zero allowed"<<std::endl;
      return 2;
    }
  }
  //check singleton in Gax: row weight = 1
  itpp::bvec bvec_zero = itpp::zeros_b(na);
  for ( int i = 0; i < Gax_row; i++){
    if ( itpp::BERC::count_errors(bvec_zero, Gax.get_row(i)) == 1){
      //      std::cout<<".";
      return 2;
    }
  }
 
  itpp::GF2mat H = common::nullSpace(Gax);
  if (debug) std::cout<<"nullSpace: H"<<H<<std::endl;
 
  //check singleton in H: row weight = 1 
  for ( int i = 0; i < na - Gax_row; i++){
    if ( itpp::BERC::count_errors(bvec_zero, H.get_row(i)) == 1){
      //      std::cout<<"+";
      return 2;
    }
  }
 
 
 
 
  //check id_Gaz
 
 
  //  if (debug) std::cout<<"rows_to_remove: "<<rows_to_remove<<std::endl;
  //  remove_rows(&H, rows_to_remove );
  if (debug) std::cout<<"alpha_Gaz"<<alpha_Gaz<<std::endl;
  Gaz=alpha_Gaz*H;
    //itpp::GF2mat(itpp::dec2bin(), false);
 
 
  //check singleton in Gaz: col weight = 0
  //  itpp::GF2mat Haz = nullSpace(Gaz);
  itpp::bvec bvec_zero_col=itpp::zeros_b(Gaz.rows());
  for ( int i = 0; i < Gaz.cols(); i++){
    if ( itpp::BERC::count_errors(bvec_zero_col, Gaz.get_col(i)) == 0){
      //      std::cout<<"+";
      return 2;
    }
  }
  
  //  Gaz = alpha_Gaz*H;
  if (debug) std::cout<<"Gaz"<<Gaz<<std::endl;
  return 0;
}
 
 
int generate_code(CSSCode & code, int debug){
  return generate_code(code.Gx, code.Gz, code.n, code.Gx_row, code.id_Gx, code.Gz_row, code.id_Gz, debug);
  }
 
 
 
 
// a version include both reduce and concatenation
// mode=1 for reduce/subsystem product
// mode=2 for concatenation
//only dz is checked cause dx is known to be tight
int product(itpp::GF2mat Gax, itpp::GF2mat Gaz, itpp::GF2mat Gbx, itpp::GF2mat Gbz,int ddax,int ddaz,int ddbx,int ddbz, int debug, int mode){
  //construct code C and calculate the distance; Compare it with the input (estimated) value
  int na=Gax.cols(),nb=Gbx.cols();//,nc=na*nb;//size of the codes
 
  itpp::GF2mat Cax = common::getC(Gax,Gaz),  Cbx = common::getC(Gbx,Gbz);//This line doesn't allow C to be empty
  itpp::GF2mat Caz = common::getC(Gax,Gaz,1),Cbz = common::getC(Gbx,Gbz,1);//This line doesn't allow C to be empty  
  
  //  Gcz=make_it_full_rank(Gcz);//not needed for calculating distance
  itpp::GF2mat Gcx,Gcz;
 
  switch ( mode ){
  case 0://reduce/subsystem product, x distance
    
    Gcx = common::kron(Gax,itpp::gf2dense_eye(nb)).concatenate_vertical(common::kron(itpp::gf2dense_eye(na),Gbx));
    Gcz=common::kron(Gaz,Gbz).concatenate_vertical(
                       common::kron(Caz,Gbz)
                       .concatenate_vertical(common::kron(Gaz,Cbz))
                       );
    break;
  case 1://reduce/subsystem product, z distance
    Gcz = common::kron(Gaz,itpp::gf2dense_eye(nb)).concatenate_vertical(common::kron(itpp::gf2dense_eye(na),Gbz));
    Gcx=common::kron(Gax,Gbx)
      .concatenate_vertical(
                common::kron(Cax,Gbx)
                .concatenate_vertical(common::kron(Gax,Cbx))
                );
    break;
  case 2://concatenation
    Gcz = common::kron(Gaz,Cbz).concatenate_vertical(common::kron(itpp::gf2dense_eye(na),Gbz));
    Gcx=common::kron(itpp::gf2dense_eye(na),Gbx).concatenate_vertical( common::kron(Gax,Cbx)   );
    break;
  case 3:
    // chain complex to two CSS codes.
  case 4:
    {
    // chain complex to two CSS codes.
    Gcx=common::kron(Gaz.transpose(), itpp::gf2dense_eye(Gbx.rows()))
      .concatenate_horizontal(common::kron(itpp::gf2dense_eye(Gax.cols()),Gbx))
      .concatenate_horizontal(common::kron(itpp::GF2mat(Gax.cols(),Gax.rows()),itpp::GF2mat(Gbx.rows(),Gbz.rows())));
    Gcx = Gcx
      .concatenate_vertical(
                common::kron(itpp::GF2mat(Gax.rows(),Gaz.rows()),itpp::GF2mat(Gbx.cols(),Gbx.rows()))
                .concatenate_horizontal(common::kron(Gax, itpp::gf2dense_eye(Gbx.cols())))
                .concatenate_horizontal(common::kron(itpp::gf2dense_eye(Gax.rows()),Gbz.transpose()))
                );
    Gcz=common::kron(itpp::gf2dense_eye(Gaz.rows()), Gbx.transpose())
      .concatenate_horizontal(common::kron(Gaz,itpp::gf2dense_eye(Gbx.cols())))
      .concatenate_horizontal(common::kron(itpp::GF2mat(Gaz.rows(),Gax.rows()),itpp::GF2mat(Gbz.cols(),Gbz.rows())));
    Gcz=Gcz
      .concatenate_vertical(
                common::kron(itpp::GF2mat( Gaz.cols(),Gaz.rows() ), itpp::GF2mat( Gbz.rows(),Gbx.rows() ))
                .concatenate_horizontal(common::kron(itpp::gf2dense_eye(Gaz.cols()),Gbz))
                .concatenate_horizontal(common::kron(Gax.transpose(),itpp::gf2dense_eye(Gbz.rows())))
                );
    
    }
    break;
 
  }  
 
 
  int flag_dist_flip=1;
  switch ( mode ){
  case 0:
  case 4:
    //x distance
    flag_dist_flip=0;
    break;
  case 1:
  case 2:
  case 3:
    // z distance
    flag_dist_flip=1;
    break;
  }
  if ( debug ){
    switch ( mode ){
    case 3:
      std::cout<<"mode (3)"<<std::endl;
      break;
    case 4:
      std::cout<<"mode (4)"<<std::endl;
      break;
    }
  }
  if (debug){common::GF2matPrint(Gcx,"Gcx"); common::GF2matPrint(Gcz,"Gcz");}
 
  if ( ! (Gcx*Gcz.transpose()).is_zero() ){
    std::cout<<"concatenation_lib: not a quantum code "<<std::endl;
    throw "not a quantum code";
  }else{
    if ( debug) std::cout<<"mode ("<<mode<<") is quantum code"<<std::endl;
  }
 
  //  Gcx=make_it_full_rank(Gcx);//not sure if I need it here
  /*  int daz = quantum_dist(Gax,Gaz,ddaz,1);
  std::cout<<"daz="<<daz<<",ddaz="<<ddaz<<std::endl;
  int dbz = quantum_dist(Gbx,Gbz,ddbz,1);
  std::cout<<"dbz="<<dbz<<",ddbz="<<ddbz<<std::endl;*/
  /*  if (is_quantum_code(Gcx,Gcz) ){
    std::cout<<"C is a quantum Code."<<std::endl;
    }*/
  //  if ( debug ) std::cout<<"Gcx"<<Gcx<<"Gcz"<<Gcz<<std::endl;
 
  switch ( flag_dist_flip ){
  case 0:
    //x distance
    {
      //  if ( debug ) std::cout<<"Gcx"<<Gcx<<"Gcz"<<Gcz<<std::endl;
      int dax=ddax,dbx=ddbx;
      int dcx = common::quantum_dist(Gcx,Gcz,dax*dbx,debug,0);//donot use estimated value ddaz and ddbz
      if (debug) std::cout<<"dax,daz,dbx,dbz = "<<ddax<<","<<ddaz<<","<<ddbx<<","<<ddbz<<","<<std::endl;    
      if (dcx == dax*dbx){
    if (debug) std::cout<<"dcx = dax*dbx = "<<dcx<<std::endl;
    return 0;
      }else if(dcx == common::INF) {
    if (debug) std::cout<<"dcx = "<<dcx<<", dax = "<<dax<<", dbx = "<<dbx<<std::endl;
    return 1;
      }else{
    if (dcx > dax*dbx) std::cout<<"PSEUDO ";
    std::cout<<common::red_text("CASE:")<<" mode ("<<mode<<") dax*dbx="<<dax*dbx<<", dcx="<<dcx;
    std::cout<<". dax,daz,dbx,dbz = "<<ddax<<","<<ddaz<<","<<ddbx<<","<<ddbz<<";";    
    std::cout<<"na,nb,nc,"<<Gax.cols()<<","<<Gbx.cols()<<","<<Gcx.cols()<<";";
    std::cout<<"ka,kb="<<Cax.rows()<<","<<Cbx.rows()<<";";
    return 2;
      }
    }
 
    break;
  case 1:
    //z distance
    {
      int daz=ddaz,dbz=ddbz;
      int dcz = common::quantum_dist(Gcx,Gcz,daz*dbz,debug,1);//donot use estimated value ddaz and ddbz
      if (debug) std::cout<<"dax,daz,dbx,dbz = "<<ddax<<","<<ddaz<<","<<ddbx<<","<<ddbz<<","<<std::endl;    
      if (dcz == daz*dbz){
    if (debug) std::cout<<"dcz = daz*dbz = "<<dcz<<std::endl;
    return 0;
      }else if(dcz == common::INF) {
    if (debug) std::cout<<"dcz = "<<dcz<<", daz = "<<daz<<", dbz = "<<dbz<<std::endl;
    return 1;
      }else{
    if (dcz > daz*dbz) std::cout<<"PSEUDO ";
    std::cout<<common::red_text("CASE:")<<" mode ("<<mode<<") daz*dbz="<<daz*dbz<<", dcz="<<dcz;
    std::cout<<". dax,daz,dbx,dbz = "<<ddax<<","<<ddaz<<","<<ddbx<<","<<ddbz<<";";
    std::cout<<"na,nb,nc="<<Gax.cols()<<","<<Gbx.cols()<<","<<Gcx.cols()<<";";    
    std::cout<<"ka,kb="<<Cax.rows()<<","<<Cbx.rows()<<";";
    return 2;
      }
    }
  }
  return 0;
  
}