diff --git a/python/optics.py b/python/optics.py index 786c3c0..9dde059 100644 --- a/python/optics.py +++ b/python/optics.py @@ -359,17 +359,12 @@ class layersample: def i(self, l, c, d): i = l * 6 + d * 3 + c - 3 return i - - #create a matrix for a single plane wave specified by k and E - # d = [dx, dy] are the x and y coordinates of the normalized direction of propagation - # k0 is the free space wave number (2 pi / lambda0) - # E is the electric field vector - - def solve1(self, d, k0, E): - - #s is the plane wave direction scaled by the refractive index - s = np.array(d) * self.n[0] + # generate the linear system corresponding to this layered sample and plane wave + # s is the direction vector scaled by the refractive index + # k0 is the free-space k-vector + # E is the field vector for the incident plane wave + def generate_linsys(self, s, k0, E): #allocate space for the matrix L = len(self.n) M = np.zeros((6*(L-1), 6*(L-1)), dtype=np.complex128) @@ -516,18 +511,33 @@ class layersample: M[ei, self.i(l, 0, 1)] = B*sz1 M[ei, self.i(l, 2, 1)] = B*s[0] ei = ei + 1 + + return [M, b] + + #create a matrix for a single plane wave specified by k and E + # d = [dx, dy] are the x and y coordinates of the normalized direction of propagation + # k0 is the free space wave number (2 pi / lambda0) + # E is the electric field vector + + def solve1(self, d, k0, E): - #store the matrix and RHS vector (for debugging) - self.M = M - self.b = b + #s is the plane wave direction scaled by the refractive index + s = np.array(d) * self.n[0] + + + #store the matrix and RHS vector (for debugging) + [self.M, self.b] = self.generate_linsys(s, k0, E) + #self.M = M + #self.b = b #evaluate the linear system - P = np.linalg.solve(M, b) + P = np.linalg.solve(self.M, self.b) #save the results (also for debugging) self.P = P #store the coefficients for each layer + L = len(self.n) # calculate the number of layers self.Pt = np.zeros((3, L), np.complex128) self.Pr = np.zeros((3, L), np.complex128) for l in range(L): diff --git a/tira/Matrix.h b/tira/Matrix.h new file mode 100644 index 0000000..f154f7c --- /dev/null +++ b/tira/Matrix.h @@ -0,0 +1,165 @@ +#include +#include +#include + +#ifndef LAPACK_COMPLEX_CUSTOM +#define LAPACK_COMPLEX_CUSTOM +#define lapack_complex_float std::complex +#define lapack_complex_double std::complex +#include "lapacke.h" +#endif + +namespace tira { + + template + class Matrix { + type* m_ptr; + int m_rows, m_cols; + + void allocate(int rows, int cols) { + if (m_ptr != NULL) + free(m_ptr); + m_rows = rows; + m_cols = cols; + m_ptr = (type*)malloc(m_rows * m_cols * sizeof(type)); + } + + public: + /// Default constructor, creates an empty matrix (generally filled as a parameter) + Matrix() { + m_rows = m_cols = 0; + m_ptr = NULL; + } + + /// Contructor, creates a rows x cols matrix with undefined values + Matrix(int rows, int cols) : Matrix() { + allocate(rows, cols); + } + + /// De-allocate data and clear the matrix for garbage collection + void del() { + m_rows = 0; + m_cols = 0; + free(m_ptr); + m_ptr = NULL; + } + + void create(int rows, int cols) { + allocate(rows, cols); + } + + int rows() { return m_rows; } + int cols() { return m_cols; } + type* get_ptr() { return m_ptr; } + + void force_square() { + if (m_rows != m_cols) + throw "Matrix must be square!"; + if (m_rows == 0 || m_cols == 0) + throw "Matrix undefined!"; + } + + inline void set(int row, int col, type val) { + m_ptr[col * m_rows + row] = val; + } + + inline type get(int row, int col) { + return m_ptr[col * m_rows + row]; + } + + /// Return the number of bytes in the matrix array + size_t bytes() { return m_rows * m_cols * sizeof(type); } + + /// Create a deep copy of the current matrix and return it + Matrix copy() { + Matrix result(m_rows, m_cols); + memcpy(result.m_ptr, m_ptr, bytes()); + return result; + } + + /// Calculate the determinant of the matrix + type det() { + force_square(); //throw an exception if the matrix isn't square + Matrix tmp = copy(); //copy the current matrix to create space for LUP decomposition + Matrix idx = getrf(tmp); //perform LUP decomposition + + // multiply all elements along the diagonal + type determinant = get(0, 0); //initialize the determinant + for (int i = 1; i < m_rows; i++) //for each remaining element of the diagonal + determinant *= tmp.get(i, i); //calculate the running factor + tmp.del(); //delete memory for the temporary matrix + + // calculate the number of row swaps + int swaps = 0; //initialize the number of row swaps + for (int i = 0; i < m_rows; i++) //iterate through each element of the swap vector + if (idx.get(i, 0) != i + 1) //if the element does not equal the row number + swaps++; //a row was swapped, so increment the swap counter + if (swaps % 2 != 0) //if there was an odd number of row swaps + determinant *= (-1); //multiply the determinant by -1 + return determinant; //return the determinant + } + + /// Calculate the trace of the matrix + type trace() { + force_square(); //throw an exception if the matrix isn't square + type trace = 0; //initialize the trace summation + for (int i = 0; i < m_rows; i++) //for each element along the diagonal + trace += get(i, i); //sum the trace + return trace; //return the matrix trace result + } + + /// Output the matrix to a CSV file + void csv(std::string filename) { + std::ofstream outfile(filename.c_str()); //open the file for writing + for (int r = 0; r < m_rows; r++) { //for each row + for (int c = 0; c < m_cols; c++) { //for each column + outfile << get(r, c); //output the value to a file + if (c != m_cols - 1) outfile << ","; //if this isn't the last column, add a comma + } + if (r != m_rows - 1) outfile << std::endl; //if this isn't the last row, end the line + } + } + + /* Static functions used to generate specialized matrices */ + + /// Static function for generating an identity matrix + static Matrix I(int dim) { + Matrix result(dim, dim); //create and allocate space for the matrix + memset(result.m_ptr, 0, dim * dim * sizeof(type)); //set all elements to zero + for (int i = 0; i < dim; i++) //for each element on the diagonal + result.set(i, i, (type)1); //set that element to 1 + return result; //return the completed matrix + } + + /// Static function for generating a geometric matrix for robustness testing + static Matrix geometric_matrix(int dim) { + Matrix result(dim, dim); //create and allocate space for the matrix + for (int r = 0; r < dim; r++) //for each row + for (int c = 0; c < dim; c++) //for each column + result.set(r, c, pow((type)(2 + r), c)); //set the geometric coefficient + return result; //return the completed matrix + } + + /// Static function analytically generates a solution (right hand side) vector for a geometric matrix + static Matrix geometric_vector(int dim) { + Matrix result(dim, 1); //generate and allocate space for the vector + for (int r = 0; r < dim; r++) + result.set(r, 0, (pow((type)(2 + r), dim) - 1) / (type)(r + 1)); + return result; //return the completed vector + } + + /* Static functions used for LAPACK and BLAS */ + + /// Calls the appropriate LAPACK GETRF function for performing LUP decomposition + static Matrix getrf(Matrix M) { + Matrix idx(M.m_rows, 1); + LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr()); + return idx; + } + static Matrix getrf(Matrix M) { + Matrix idx(M.m_rows, 1); + LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr()); + return idx; + } + }; +} \ No newline at end of file -- libgit2 0.21.4