function mtx_A = coefficient_A(n1, n2, n3, sigma, vec_k)

cubic_length      = 1;
period_vec1       = cubic_length * [ 1; 0; 0 ];
period_vec2       = cubic_length * [ 0; 1; 0 ];
period_vec3       = cubic_length * [ 0; 0; 1 ];

length_x          = cubic_length;
length_y          = cubic_length;
length_z          = cubic_length;

delta              = [ length_x / n1; length_y / n2; length_z / n3];

 mtx_G            = Sparse_mtx_G(n1,n2,n3,delta,period_vec1,period_vec2,period_vec3,vec_k);
 mtx_A            = mtx_G' * mtx_G + sigma * speye(3*n1*n2*n3);
 
end

% ===========================

 function [mtx_G] = Sparse_mtx_G(n_1,n_2,n_3,delta,period_vec1,period_vec2,period_vec3,vec_k)

     n12  = n_1 * n_2;
     n123 = n12 * n_3;
     
     [i_row_1, j_col_1, val_1] = Gamma_1(n_1, n_2, n_3, period_vec1, vec_k);
     [i_row_2, j_col_2, val_2] = Gamma_2(n_1, n_2, n_3, period_vec2, vec_k);
     [i_row_3, j_col_3, val_3] = Gamma_3(n_1,n_2,n_3, period_vec3, vec_k);
     
     leng_1 = length(i_row_1);
     leng_2 = length(i_row_2);
     leng_3 = length(i_row_3);
     leng   = 2 * ( leng_1  + leng_2 + leng_3 );
     i_row  = zeros(leng,1);
     j_col  = zeros(leng,1);
     val    = zeros(leng,1);
     
     val_1  = val_1 / delta(1);
     val_2  = val_2 / delta(2);
     val_3  = val_3 / delta(3);
     
%     mtx_C3 = sparse(i_row_3, j_col_3, val_3);
     
%    First block rows     
     i_row(1:leng_3,1) = i_row_3;
     j_col(1:leng_3,1) = n123 + j_col_3;
     val  (1:leng_3,1) = val_3;
     
     ii                   = leng_3 + leng_2;
     i_row(leng_3+1:ii,1) = i_row_2;
     j_col(leng_3+1:ii,1) = 2 * n123 + j_col_2;
     val  (leng_3+1:ii,1) = - val_2;
% Second block rows     
     i2               = ii + leng_3;
     i_row(ii+1:i2,1) = n123 + i_row_3;
     j_col(ii+1:i2,1) = j_col_3;
     val  (ii+1:i2,1) = - val_3;
     
     ii               = i2 + leng_1;
     i_row(i2+1:ii,1) = n123 + i_row_1;
     j_col(i2+1:ii,1) = 2 * n123 + j_col_1;
     val  (i2+1:ii,1) = val_1;
% Third block rows     
     i2               = ii + leng_2;
     i_row(ii+1:i2,1) = 2 * n123 + i_row_2;
     j_col(ii+1:i2,1) = j_col_2;
     val  (ii+1:i2,1) = val_2;
     
     ii               = i2 + leng_1;
     i_row(i2+1:ii,1) = 2 * n123 + i_row_1;
     j_col(i2+1:ii,1) = n123 + j_col_1;
     val  (i2+1:ii,1) = - val_1;
     
     mtx_G            = sparse(i_row, j_col, val);
     
 end
 
 % =========================================
 function [i_row, j_col, val] = Gamma_1(n_1, n_2, n_3, period_vec1, vec_k)
 
mtx_C_1 = gen_mtx_C_1(n_1,period_vec1,vec_k);

[i_row_tmp, j_col_tmp, val_tmp] = find(mtx_C_1);
leng                            = length(i_row_tmp);
i_row                           = zeros(leng*n_2*n_3,1);
j_col                           = zeros(leng*n_2*n_3,1);
val                             = zeros(leng*n_2*n_3,1);
for ii = 1:n_2*n_3
    jj                    = ( ii - 1 ) * n_1;
    kk                    = ( ii - 1 ) * leng;
    i_row(kk+1:kk+leng,1) = jj + i_row_tmp;
    j_col(kk+1:kk+leng,1) = jj + j_col_tmp;
    val  (kk+1:kk+leng,1) = val_tmp;
end

 end
 
 % ========================================
 
 function [mtx_C_1] = gen_mtx_C_1(n_1,period_vec1,vec_k)
%
% Generate matrix K_1^{*}
%
     mtx_C_1 = -eye(n_1);

     for ii = 1:n_1-1
        mtx_C_1(ii,ii+1) = 1;
     end

     tmp            = 2 * pi * (vec_k' * period_vec1) * 1i;
     mtx_C_1(n_1,1) = exp(tmp);
     
 end
 
 % ========================================
 function [i_row, j_col, val] = Gamma_2(n_1, n_2, n_3, period_vec2, vec_k)
 
[i_row_tmp, j_col_tmp, val_tmp] = gen_mtx_C_2(n_1,n_2,period_vec2,vec_k); 
leng                            = length(i_row_tmp);
i_row                           = zeros(leng*n_3,1);
j_col                           = zeros(leng*n_3,1);
val                             = zeros(leng*n_3,1);
for ii = 1:n_3
    jj                    = ( ii - 1 ) * n_1 * n_2;
    kk                    = ( ii - 1 ) * leng;
    i_row(kk+1:kk+leng,1) = jj + i_row_tmp;
    j_col(kk+1:kk+leng,1) = jj + j_col_tmp;
    val  (kk+1:kk+leng,1) = val_tmp;
end

 end
 
 % =================================================
 
 function [i_row, j_col, val] = gen_mtx_C_2(n_1,n_2,period_vec2,vec_k)
%
% Generate matrix K_2^{*}
%
     i_row = zeros(2*n_1*n_2,1);
     j_col = zeros(2*n_1*n_2,1);
     val   = zeros(2*n_1*n_2,1);
     for jj = 1:n_2-1
         ii                   = ( jj - 1 ) * n_1;
         kk                   = 2 * ii;
         i_row(kk+1:kk+n_1,1) = ii+1:ii+n_1;
         j_col(kk+1:kk+n_1,1) = ii+1:ii+n_1;
         val  (kk+1:kk+n_1,1) = -ones(n_1,1);
         
         i_row(kk+n_1+1:kk+2*n_1,1) = ii+1:ii+n_1;
         j_col(kk+n_1+1:kk+2*n_1,1) = ii+n_1+1:ii+2*n_1;
         val  (kk+n_1+1:kk+2*n_1,1) = ones(n_1,1);
     end
   
     ii                   = ( n_2 - 1 ) * n_1;
     kk                   = 2 * ii;
     i_row(kk+1:kk+n_1,1) = ii+1:ii+n_1;
     j_col(kk+1:kk+n_1,1) = ii+1:ii+n_1;
     val  (kk+1:kk+n_1,1) = -ones(n_1,1);
         
     i_row(kk+n_1+1:kk+2*n_1,1) = ii+1:ii+n_1;
     j_col(kk+n_1+1:kk+2*n_1,1) = 1:n_1;
     tmp                        = 2 * pi * (vec_k' * period_vec2) * 1i;
     val  (kk+n_1+1:kk+2*n_1,1) = exp(tmp) * ones(n_1,1);
         
 end
 
 % =================================================
 
 function [i_row, j_col, val] = Gamma_3(n_1,n_2,n_3, period_vec3, vec_k)
 
     n12   = n_1 * n_2; 
     i_row = zeros(2*n_1*n_2*n_3,1);
     j_col = zeros(2*n_1*n_2*n_3,1);
     val   = zeros(2*n_1*n_2*n_3,1);
     for jj = 1:n_3-1
         ii                   = ( jj - 1 ) * n12;
         kk                   = 2 * ii;
         i_row(kk+1:kk+n12,1) = ii+1:ii+n12;
         j_col(kk+1:kk+n12,1) = ii+1:ii+n12;
         val  (kk+1:kk+n12,1) = -ones(n12,1);
         
         i_row(kk+n12+1:kk+2*n12,1) = ii+1:ii+n12;
         j_col(kk+n12+1:kk+2*n12,1) = ii+n12+1:ii+2*n12;
         val  (kk+n12+1:kk+2*n12,1) = ones(n12,1);
     end
   
     ii                   = ( n_3 - 1 ) * n12;
     kk                   = 2 * ii;
     i_row(kk+1:kk+n12,1) = ii+1:ii+n12;
     j_col(kk+1:kk+n12,1) = ii+1:ii+n12;
     val  (kk+1:kk+n12,1) = -ones(n12,1);
         
     i_row(kk+n12+1:kk+2*n12,1) = ii+1:ii+n12;
     j_col(kk+n12+1:kk+2*n12,1) = 1:n12;
     tmp                        = 2 * pi * (vec_k' * period_vec3) * 1i;
     val  (kk+n12+1:kk+2*n12,1) = exp(tmp) * ones(n12,1);    
     
 end
 
 % ===============================================
 