# Path: sucemfem/Utilities/LinalgSolvers.py
# def solve_sparse_system ( A, b, preconditioner_type='ilu' ):
#     """
#     This function solves the sparse linear system Ax = b for A a scipy sparse matrix, and b a numpy or scipy array
#     
#     By default an incomplete LU preconditioner is used with the bicgstab iterative solver
#     
#     @param A: a square matrix 
#     @param b: the RHS vector
#     @param preconditioner_type: Preconditioner type string
#     """
#                     
#     solver = BiCGStabSolver ( A, preconditioner_type )
# #    solver = BiCGStabSolver ( A ) #removed preconditioner -> 2D Waveguide works!    
#     x = solver.solve(b)
#     return x
# 
# Path: sucemfem/Utilities/Converters.py
# def dolfin_ublassparse_to_scipy_csr ( A, dtype=None, imagify=False ):
#     """
#     convert a DOLFIN uBLASSparseMatrix to a scipy.sparse.csr_matrix()
#     
#     @param A: a DOLFIN uBLASSparseMatrix
#     @param dtype: the numpy data type to use to store the matrix
#     @param imagify: multiply the original matrix data by 1j
#     """
#     import scipy.sparse
#     # get the sparse data from the input matrix
#     (row,col,data) = A.data()   # get sparse data
#     col = np.intc(col)
#     row = np.intc(row)
#     n = A.size(0)
#     if imagify: data = data*1j
#     A_sp = scipy.sparse.csr_matrix( (data,col,row), shape=(n,n), dtype=dtype)
#     
#     return A_sp


import numpy as N
import pylab as P
import scipy.sparse
import csv
import sys
from dolfin import *
from sucemfem.Utilities.LinalgSolvers import solve_sparse_system
from sucemfem.Utilities.Converters import dolfin_ublassparse_to_scipy_csr
    print "Setting matrix entries"
    t = Timer("Set matrix entries")
    Beta = Beta_km(f)
    ko = k_o(f)
    coefficient_H_C = 1j*Beta*2/N.sqrt(a*b)*N.sqrt(ko*Z_0/Beta)
    coefficient_A_E = N.sqrt(a/b)*N.sqrt(ko*Z_0/Beta)
    
    for k in range(num_ports):
        for m in range(M):
            j = k*M + m
            
            B_km = Beta_km(f, m+1)
            AD[j,j] = -N.sqrt(a/b)*N.sqrt(ko*Z_0/B_km)
            coefficients_C[j] = 1j*B_km*2/N.sqrt(a*b)*N.sqrt(ko*Z_0/B_km)
            
    
    F_sp = dolfin_ublassparse_to_scipy_csr ( F ) # scipy.sparse.csr_matrix ( F.array() )
    C_sp = scipy.sparse.csr_matrix ( N.dot(C, N.diag(coefficients_C)) )

#    LHS[num_ports*M:,num_ports*M:] = F.array()[:,:]
#    LHS[num_ports*M:,:num_ports*M] = N.dot(C, N.diag(coefficients_C))
    LHS = scipy.sparse.vstack( ( AD, scipy.sparse.hstack( (C_sp, F_sp ), 'csr') ), 'csr' )
    
    RHS[0,0] = -LHS[0,0]
    RHS[num_ports*M:,0] = coefficient_H_C*H_inc.array()[:]
    t.stop()
    
    print "Solving (sparse)"
    t = Timer("Sparse Solve")
    

# Path: sucemfem/Utilities/LinalgSolvers.py
# def solve_sparse_system ( A, b, preconditioner_type='ilu' ):
#     """
#     This function solves the sparse linear system Ax = b for A a scipy sparse matrix, and b a numpy or scipy array
#     
#     By default an incomplete LU preconditioner is used with the bicgstab iterative solver
#     
#     @param A: a square matrix 
#     @param b: the RHS vector
#     @param preconditioner_type: Preconditioner type string
#     """
#                     
#     solver = BiCGStabSolver ( A, preconditioner_type )
# #    solver = BiCGStabSolver ( A ) #removed preconditioner -> 2D Waveguide works!    
#     x = solver.solve(b)
#     return x
# 
# Path: sucemfem/Utilities/Converters.py
# def dolfin_ublassparse_to_scipy_csr ( A, dtype=None, imagify=False ):
#     """
#     convert a DOLFIN uBLASSparseMatrix to a scipy.sparse.csr_matrix()
#     
#     @param A: a DOLFIN uBLASSparseMatrix
#     @param dtype: the numpy data type to use to store the matrix
#     @param imagify: multiply the original matrix data by 1j
#     """
#     import scipy.sparse
#     # get the sparse data from the input matrix
#     (row,col,data) = A.data()   # get sparse data
#     col = np.intc(col)
#     row = np.intc(row)
#     n = A.size(0)
#     if imagify: data = data*1j
#     A_sp = scipy.sparse.csr_matrix( (data,col,row), shape=(n,n), dtype=dtype)
#     
#     return A_sp


# Path: examples/2D_waveguide/h_plane_hollow_rectangular_S_parameters.py
import numpy as N
import pylab as P
import scipy.sparse
import csv
import sys
from dolfin import *
from sucemfem.Utilities.LinalgSolvers import solve_sparse_system
from sucemfem.Utilities.Converters import dolfin_ublassparse_to_scipy_csr

    print "Setting matrix entries"
    t = Timer("Set matrix entries")
    Beta = Beta_km(f)
    ko = k_o(f)
    coefficient_H_C = 1j*Beta*2/N.sqrt(a*b)*N.sqrt(ko*Z_0/Beta)
    coefficient_A_E = N.sqrt(a/b)*N.sqrt(ko*Z_0/Beta)
    
    for k in range(num_ports):
        for m in range(M):
            j = k*M + m
            
            B_km = Beta_km(f, m+1)
            AD[j,j] = -N.sqrt(a/b)*N.sqrt(ko*Z_0/B_km)
            coefficients_C[j] = 1j*B_km*2/N.sqrt(a*b)*N.sqrt(ko*Z_0/B_km)
            
    
    F_sp = dolfin_ublassparse_to_scipy_csr ( F ) # scipy.sparse.csr_matrix ( F.array() )
    C_sp = scipy.sparse.csr_matrix ( N.dot(C, N.diag(coefficients_C)) )

#    LHS[num_ports*M:,num_ports*M:] = F.array()[:,:]
#    LHS[num_ports*M:,:num_ports*M] = N.dot(C, N.diag(coefficients_C))
    LHS = scipy.sparse.vstack( ( AD, scipy.sparse.hstack( (C_sp, F_sp ), 'csr') ), 'csr' )
    
    RHS[0,0] = -LHS[0,0]
    RHS[num_ports*M:,0] = coefficient_H_C*H_inc.array()[:]
    t.stop()
    
    print "Solving (sparse)"
    t = Timer("Sparse Solve")
    

    BE = solve_sparse_system (LHS, RHS)