#! /usr/bin/env python

# -*- coding: utf-8 -*-
"""
Created on Tue Sep 30 11:25:36 2014

@author: hans-werner
"""

def hw3q3():
    """ Investigate the semi-positive definiteness of system matrices """
        
    import numpy as np
    import scipy.linalg as la
    
    #
    # Define the mesh
    #
    a = 0.0
    b = 1.0
    N = 5
    x = np.linspace(a,b,N+1)
    
    #
    #  Set a 3 point quadrature rule on the reference interval [0,1].
    #
    ng = 3

    xg = np.array ( ( \
    0.112701665379258311482073460022, \
    0.5, \
    0.887298334620741688517926539978 ) )

    wg = np.array ( ( \
    5.0 / 18.0, \
    8.0 / 18.0, \
    5.0 / 18.0 ) )
    
    
    #
    # Compute the three system matrices 
    #
    K = np.zeros([N+1,N+1])
    L = np.zeros([N+1,N+1])
    M = np.zeros([N+1,N+1])

    #    
    # Loop over the elements
    #
    for e in range(N):
        
        xl = x[e]
        xr = x[e+1]
        
        #
        #  Consider quadrature point Q: (0, 1, 2 ) in element E.
        #
        for q in range ( 0, ng ):
            #
            #  Map XG and WG from [0,1] to
            #      XQ and QQ in [XL,XR].
            #
            xq = xl + xg[q] * ( xr - xl )
            wq = wg[q] * ( xr - xl )
        
            #
            #  Consider the I-th test function PHI(I,X) and 
            #  its derivative PHI'(I,X).
            #
            for i_local in range ( 0, 2 ):
                i = i_local + e

                if ( i_local == 0 ):
                    phii = ( xq - xr ) / ( xl - xr )
                    phiip = 1.0 / ( xl - xr )
                else:
                    phii = ( xq - xl ) / ( xr - xl )
                    phiip = 1.0 / ( xr - xl )

        
                #
                #  Consider the J-th basis function PHI(J,X) and 
                #  its derivative PHI'(J,X). 
                #
                for j_local in range ( 0, 2 ):
                    j = j_local + e
            
                    if ( j_local == 0 ):
                        phij = ( xq - xr ) / ( xl - xr )  
                        phijp = 1.0 / ( xl - xr )
                    else:
                        phij = ( xq - xl ) / ( xr - xl )  
                        phijp = 1.0 / ( xr - xl )
            
                    # 
                    # Update the system matrices
                    # 
                    K[i][j] = K[i][j] + wq*phiip*phijp
                    L[i][j] = L[i][j] + wq*phii*phijp
                    M[i][j] = M[i][j] + wq*phii*phij
    # 
    # Combine the matrices to form the system matrix
    A = K + L + M
    
    #
    # Compute the matrices' eigenvalues
    #
    lmd_K,_ = la.eigh(K)
    lmd_L,_ = la.eig(L)
    lmd_M,_ = la.eigh(M)
    lmd_A,_ = la.eig(A)
    # 
    # Print the Matrices and their eigenvalues
    #  
    print " K : Matrix\n"
    for i in range(N+1):
        for j in range(N+1):
            print "%3.0f" % (K[i][j]),
        print "\n"         
            
    print "K: Eigenvalues\n"
    for i in range(N+1):    
        print "%2.3f" % (lmd_K[i]),
    print "\n\n"
    
    print "L : Matrix \n"
    for i in range(N+1):
        for j in range(N+1):
            print "%4.1f" % (L[i][j]),
        print "\n"

    print "L: Eigenvalues\n"
    for i in range(N+1):    
        print "%1.0f + %1.2fi, " % (lmd_L[i].real,lmd_L[i].imag),
    print "\n\n"
    
    print "M : Matrix\n"
    for i in range(N+1):
        for j in range(N+1):
            print "%4.5f" % (M[i][j]),
        print "\n" 
    
    print "M: Eigenvalues\n"
    for i in range(N+1):    
        print "%2.3f" % (lmd_M[i]),
    print "\n\n"
    
    print "A : Matrix\n"
    for i in range(N+1):
        for j in range(N+1):
            print "%5.3f" % (A[i][j]),
        print "\n" 
    
    print "A: Eigenvalues\n"
    for i in range(N+1):    
        print "%2.3f + %1.0f i," % (lmd_A[i].real,lmd_A[i].imag),
    print "\n\n"    
#    
#  If this script is called directly, then run it as a program.
#
if ( __name__ == '__main__' ):
  hw3q3( )