

function [ GtbsTri, CtbsTri, BtbsTri ] = tri ( Gg, Cg, Bg, Vtg, GgSBB_tot, CgSBB_tot, GXconn, CXconn , SH ) 

global s0 p0 layer maxLevel numBlocks max_row_idx max_clm_idx q0 
global ugH ucH  ugV ucV  ugHX ucHX  ugVX ucVX  ugx ucx

%%%%%%%%%%%%%%%%  TRIGONALIZE %%%%%%%%%%%%%%%%%%%%%%%%%%%

 BtbsTri.xq = sparse ( 2*q0*layer*numBlocks*(maxLevel+1) , p0*numBlocks ) ; 
 idx=1 ;
 for i = 1 : numBlocks 
     BtbsTri.xq( (i-1)*(maxLevel+1)*q0*layer + 1 : (i-1)*(maxLevel+1)*q0*layer + q0*layer , (idx-1)*p0+1 : (idx-1)*p0+p0 ) = Bg(idx).xq ;
     idx = idx + 1 ;
 end
 
 
 idx=1;
 for i = numBlocks+1 : 2*numBlocks 
     BtbsTri.xq( (i-1)*(maxLevel+1)*q0*layer + 1 : (i-1)*(maxLevel+1)*q0*layer + q0*layer , (idx-1)*p0+1 : (idx-1)*p0 + p0 ) = Bg(idx).xq ;
     idx = idx + 1 ;
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for i=1:numBlocks
      for j=1:maxLevel
            GgSBB(i,j).xq = Vtg(i).x'*GgSBB_tot(j).x*Vtg(i).x ;
            CgSBB(i,j).xq = Vtg(i).x'*CgSBB_tot(j).x*Vtg(i).x ;
        end
 end
 idx=1;     
 for i=1:numBlocks
	for j=1:maxLevel
     	tmp1 ( i , j ).xq = GgSBB(i,j).xq * SH ( idx, 1 ) ;
        tmp2 ( i , j ).xq = CgSBB(i,j).xq * SH ( idx ,1 );
        idx = idx + 1 ;
	end
 end
    

%%%%%%%%%%%%%%%  Begin Iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 for i=1:numBlocks
         G01(i).xq = sparse ( q0*layer*(maxLevel+1) , q0*layer*(maxLevel+1) ) ;
         C01(i).xq = sparse ( q0*layer*(maxLevel+1) , q0*layer*(maxLevel+1) ) ;
 end


 for i=1:numBlocks
	for j=1:maxLevel+1
     		G01(i).xq (  (j-1)*layer*q0+1 :(j-1)*layer*q0+layer*q0 ,  (j-1)*layer*q0+1 : (j-1)*layer*q0 + layer*q0 ) = Gg(i).xq ;
     		C01(i).xq (  (j-1)*layer*q0+1 :(j-1)*layer*q0+layer*q0 ,  (j-1)*layer*q0+1 : (j-1)*layer*q0 + layer*q0 ) = Cg(i).xq ;
	end
  end


  for i=1:numBlocks
	for j=1:maxLevel
      		G01(i).xq ( q0 * layer + (j-1)*q0*layer + 1 : q0 * layer + (j-1)*q0*layer + q0*layer   , 1 : q0*layer  )= tmp1(i,j).xq ; 
      		C01(i).xq ( q0 * layer + (j-1)*q0*layer + 1 : q0 * layer + (j-1)*q0*layer + q0*layer   , 1 : q0*layer  )= tmp2(i,j).xq ; 
	end
  end


  %%%%%%%%%%%%% Triangularize matrix

  Gtbs_ul.xq = sparse ( (maxLevel+1)*numBlocks*layer*q0 , (maxLevel+1)*numBlocks*layer*q0 ) ;
  Gtbs_ur.xq = sparse ( (maxLevel+1)*numBlocks*layer*q0 , (maxLevel+1)*numBlocks*layer*q0 ) ;
  Gtbs_lr.xq = sparse ( (maxLevel+1)*numBlocks*layer*q0 , (maxLevel+1)*numBlocks*layer*q0 ) ;

  Ctbs_ul.xq = sparse ( (maxLevel+1)*numBlocks*layer*q0 , (maxLevel+1)*numBlocks*layer*q0 ) ;
  Ctbs_ur.xq = sparse ( (maxLevel+1)*numBlocks*layer*q0 , (maxLevel+1)*numBlocks*layer*q0 ) ;
  Ctbs_lr.xq = sparse ( (maxLevel+1)*numBlocks*layer*q0 , (maxLevel+1)*numBlocks*layer*q0 ) ;


  for k=1:numBlocks
      for l=k:numBlocks
          if ( k==l )
              Gtbs_ul.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (l-1)*layer*q0*(maxLevel+1) + 1 : (l-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = G01(k).xq ;
              Ctbs_ul.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (l-1)*layer*q0*(maxLevel+1) + 1 : (l-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = C01(k).xq ;

          else 
              Gtbs_ul.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (l-1)*layer*q0*(maxLevel+1) + 1 : (l-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = GXconn(k,l).xq ;
              Ctbs_ul.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (l-1)*layer*q0*(maxLevel+1) + 1 : (l-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = CXconn(k,l).xq ;
          end
      end
  end

  for k=1:numBlocks
      for l=1:k-1
             Gtbs_ur.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (l-1)*layer*q0*(maxLevel+1) + 1 : (l-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = GXconn(k,l).xq ;
             Ctbs_ur.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (l-1)*layer*q0*(maxLevel+1) + 1 : (l-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = CXconn(k,l).xq ;
      end
 end

 for k=1:numBlocks
      for m=1:numBlocks
        if ( k==m )
              Gtbs_lr.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (m-1)*layer*q0*(maxLevel+1) + 1 : (m-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = G01(k).xq ;
              Ctbs_lr.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (m-1)*layer*q0*(maxLevel+1) + 1 : (m-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = C01(k).xq ;
           else 
              Gtbs_lr.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (m-1)*layer*q0*(maxLevel+1) + 1 : (m-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = GXconn(k,m).xq ;
             Ctbs_lr.xq ( (k-1)*layer*q0*(maxLevel+1) + 1 : (k-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) , (m-1)*layer*q0*(maxLevel+1) + 1 : (m-1)*layer*q0*(maxLevel+1) + layer*q0*(maxLevel+1) ) = CXconn(k,m).xq ;
           end 
       end
 end

 
 GtbsTri.xq = sparse ( 2 * (maxLevel+1) * numBlocks * layer * q0 , 2 * (maxLevel+1) * numBlocks * layer * q0 ) ;
 CtbsTri.xq = sparse ( 2 * (maxLevel+1) * numBlocks * layer * q0 , 2 * (maxLevel+1) * numBlocks * layer * q0 ) ;

 GtbsTri.xq ( 1 : q0*layer*(maxLevel+1)*numBlocks , 1 : q0*layer*(maxLevel+1)*numBlocks ) = Gtbs_ul.xq ;
 GtbsTri.xq ( 1 : q0*layer*(maxLevel+1)*numBlocks , q0*layer*(maxLevel+1)*numBlocks + 1 : 2*q0*layer*(maxLevel+1)*numBlocks ) = Gtbs_ur.xq ;
 GtbsTri.xq ( q0*layer*(maxLevel+1)*numBlocks + 1 : 2*q0*layer*(maxLevel+1)*numBlocks , q0*layer*(maxLevel+1)*numBlocks + 1 : 2*q0*layer*(maxLevel+1)*numBlocks ) = Gtbs_lr.xq;

 CtbsTri.xq ( 1 : q0*layer*(maxLevel+1)*numBlocks , 1 : q0*layer*(maxLevel+1)*numBlocks ) = Ctbs_ul.xq ;
 CtbsTri.xq ( 1 : q0*layer*(maxLevel+1)*numBlocks , q0*layer*(maxLevel+1)*numBlocks + 1 : 2*q0*layer*(maxLevel+1)*numBlocks ) = Ctbs_ur.xq ;
 CtbsTri.xq ( q0*layer*(maxLevel+1)*numBlocks + 1 : 2*q0*layer*(maxLevel+1)*numBlocks , q0*layer*(maxLevel+1)*numBlocks + 1 : 2*q0*layer*(maxLevel+1)*numBlocks ) = Ctbs_lr.xq;
    

