codebase / stimlib

Commit 9333bc22b19f525f7d197d7c4a4dc48d0cd3b3b1

Authored by David Mayerich
1 parent 057bed5a

started TIRA library

Showing 2 changed files with 189 additions and 14 deletions   Show diff stats
Inline Side-by-side
python/optics.py
  Show/Hide comments   View file @9333bc2
... ... @@ -359,17 +359,12 @@ class layersample:
359 359 def i(self, l, c, d):
360 360 i = l * 6 + d * 3 + c - 3
361 361 return i
362   -
363   - #create a matrix for a single plane wave specified by k and E
364   - # d = [dx, dy] are the x and y coordinates of the normalized direction of propagation
365   - # k0 is the free space wave number (2 pi / lambda0)
366   - # E is the electric field vector
367   -
368   - def solve1(self, d, k0, E):
369   -
370   - #s is the plane wave direction scaled by the refractive index
371   - s = np.array(d) * self.n[0]
372 362  
  363 + # generate the linear system corresponding to this layered sample and plane wave
  364 + # s is the direction vector scaled by the refractive index
  365 + # k0 is the free-space k-vector
  366 + # E is the field vector for the incident plane wave
  367 + def generate_linsys(self, s, k0, E):
373 368 #allocate space for the matrix
374 369 L = len(self.n)
375 370 M = np.zeros((6*(L-1), 6*(L-1)), dtype=np.complex128)
... ... @@ -516,18 +511,33 @@ class layersample:
516 511 M[ei, self.i(l, 0, 1)] = B*sz1
517 512 M[ei, self.i(l, 2, 1)] = B*s[0]
518 513 ei = ei + 1
  514 +
  515 + return [M, b]
  516 +
  517 + #create a matrix for a single plane wave specified by k and E
  518 + # d = [dx, dy] are the x and y coordinates of the normalized direction of propagation
  519 + # k0 is the free space wave number (2 pi / lambda0)
  520 + # E is the electric field vector
  521 +
  522 + def solve1(self, d, k0, E):
519 523  
520   - #store the matrix and RHS vector (for debugging)
521   - self.M = M
522   - self.b = b
  524 + #s is the plane wave direction scaled by the refractive index
  525 + s = np.array(d) * self.n[0]
  526 +
  527 +
  528 + #store the matrix and RHS vector (for debugging)
  529 + [self.M, self.b] = self.generate_linsys(s, k0, E)
  530 + #self.M = M
  531 + #self.b = b
523 532  
524 533 #evaluate the linear system
525   - P = np.linalg.solve(M, b)
  534 + P = np.linalg.solve(self.M, self.b)
526 535  
527 536 #save the results (also for debugging)
528 537 self.P = P
529 538  
530 539 #store the coefficients for each layer
  540 + L = len(self.n) # calculate the number of layers
531 541 self.Pt = np.zeros((3, L), np.complex128)
532 542 self.Pr = np.zeros((3, L), np.complex128)
533 543 for l in range(L):
... ...
tira/Matrix.h 0 → 100644
  Show/Hide comments   View file @9333bc2
  1 +#include <iostream>
  2 +#include <fstream>
  3 +#include <complex>
  4 +
  5 +#ifndef LAPACK_COMPLEX_CUSTOM
  6 +#define LAPACK_COMPLEX_CUSTOM
  7 +#define lapack_complex_float std::complex<float>
  8 +#define lapack_complex_double std::complex<double>
  9 +#include "lapacke.h"
  10 +#endif
  11 +
  12 +namespace tira {
  13 +
  14 + template <typename type>
  15 + class Matrix {
  16 + type* m_ptr;
  17 + int m_rows, m_cols;
  18 +
  19 + void allocate(int rows, int cols) {
  20 + if (m_ptr != NULL)
  21 + free(m_ptr);
  22 + m_rows = rows;
  23 + m_cols = cols;
  24 + m_ptr = (type*)malloc(m_rows * m_cols * sizeof(type));
  25 + }
  26 +
  27 + public:
  28 + /// Default constructor, creates an empty matrix (generally filled as a parameter)
  29 + Matrix() {
  30 + m_rows = m_cols = 0;
  31 + m_ptr = NULL;
  32 + }
  33 +
  34 + /// Contructor, creates a rows x cols matrix with undefined values
  35 + Matrix(int rows, int cols) : Matrix() {
  36 + allocate(rows, cols);
  37 + }
  38 +
  39 + /// De-allocate data and clear the matrix for garbage collection
  40 + void del() {
  41 + m_rows = 0;
  42 + m_cols = 0;
  43 + free(m_ptr);
  44 + m_ptr = NULL;
  45 + }
  46 +
  47 + void create(int rows, int cols) {
  48 + allocate(rows, cols);
  49 + }
  50 +
  51 + int rows() { return m_rows; }
  52 + int cols() { return m_cols; }
  53 + type* get_ptr() { return m_ptr; }
  54 +
  55 + void force_square() {
  56 + if (m_rows != m_cols)
  57 + throw "Matrix must be square!";
  58 + if (m_rows == 0 || m_cols == 0)
  59 + throw "Matrix undefined!";
  60 + }
  61 +
  62 + inline void set(int row, int col, type val) {
  63 + m_ptr[col * m_rows + row] = val;
  64 + }
  65 +
  66 + inline type get(int row, int col) {
  67 + return m_ptr[col * m_rows + row];
  68 + }
  69 +
  70 + /// Return the number of bytes in the matrix array
  71 + size_t bytes() { return m_rows * m_cols * sizeof(type); }
  72 +
  73 + /// Create a deep copy of the current matrix and return it
  74 + Matrix copy() {
  75 + Matrix result(m_rows, m_cols);
  76 + memcpy(result.m_ptr, m_ptr, bytes());
  77 + return result;
  78 + }
  79 +
  80 + /// Calculate the determinant of the matrix
  81 + type det() {
  82 + force_square(); //throw an exception if the matrix isn't square
  83 + Matrix tmp = copy(); //copy the current matrix to create space for LUP decomposition
  84 + Matrix<int> idx = getrf(tmp); //perform LUP decomposition
  85 +
  86 + // multiply all elements along the diagonal
  87 + type determinant = get(0, 0); //initialize the determinant
  88 + for (int i = 1; i < m_rows; i++) //for each remaining element of the diagonal
  89 + determinant *= tmp.get(i, i); //calculate the running factor
  90 + tmp.del(); //delete memory for the temporary matrix
  91 +
  92 + // calculate the number of row swaps
  93 + int swaps = 0; //initialize the number of row swaps
  94 + for (int i = 0; i < m_rows; i++) //iterate through each element of the swap vector
  95 + if (idx.get(i, 0) != i + 1) //if the element does not equal the row number
  96 + swaps++; //a row was swapped, so increment the swap counter
  97 + if (swaps % 2 != 0) //if there was an odd number of row swaps
  98 + determinant *= (-1); //multiply the determinant by -1
  99 + return determinant; //return the determinant
  100 + }
  101 +
  102 + /// Calculate the trace of the matrix
  103 + type trace() {
  104 + force_square(); //throw an exception if the matrix isn't square
  105 + type trace = 0; //initialize the trace summation
  106 + for (int i = 0; i < m_rows; i++) //for each element along the diagonal
  107 + trace += get(i, i); //sum the trace
  108 + return trace; //return the matrix trace result
  109 + }
  110 +
  111 + /// Output the matrix to a CSV file
  112 + void csv(std::string filename) {
  113 + std::ofstream outfile(filename.c_str()); //open the file for writing
  114 + for (int r = 0; r < m_rows; r++) { //for each row
  115 + for (int c = 0; c < m_cols; c++) { //for each column
  116 + outfile << get(r, c); //output the value to a file
  117 + if (c != m_cols - 1) outfile << ","; //if this isn't the last column, add a comma
  118 + }
  119 + if (r != m_rows - 1) outfile << std::endl; //if this isn't the last row, end the line
  120 + }
  121 + }
  122 +
  123 + /* Static functions used to generate specialized matrices */
  124 +
  125 + /// Static function for generating an identity matrix
  126 + static Matrix I(int dim) {
  127 + Matrix result(dim, dim); //create and allocate space for the matrix
  128 + memset(result.m_ptr, 0, dim * dim * sizeof(type)); //set all elements to zero
  129 + for (int i = 0; i < dim; i++) //for each element on the diagonal
  130 + result.set(i, i, (type)1); //set that element to 1
  131 + return result; //return the completed matrix
  132 + }
  133 +
  134 + /// Static function for generating a geometric matrix for robustness testing
  135 + static Matrix geometric_matrix(int dim) {
  136 + Matrix result(dim, dim); //create and allocate space for the matrix
  137 + for (int r = 0; r < dim; r++) //for each row
  138 + for (int c = 0; c < dim; c++) //for each column
  139 + result.set(r, c, pow((type)(2 + r), c)); //set the geometric coefficient
  140 + return result; //return the completed matrix
  141 + }
  142 +
  143 + /// Static function analytically generates a solution (right hand side) vector for a geometric matrix
  144 + static Matrix geometric_vector(int dim) {
  145 + Matrix result(dim, 1); //generate and allocate space for the vector
  146 + for (int r = 0; r < dim; r++)
  147 + result.set(r, 0, (pow((type)(2 + r), dim) - 1) / (type)(r + 1));
  148 + return result; //return the completed vector
  149 + }
  150 +
  151 + /* Static functions used for LAPACK and BLAS */
  152 +
  153 + /// Calls the appropriate LAPACK GETRF function for performing LUP decomposition
  154 + static Matrix<int> getrf(Matrix<float> M) {
  155 + Matrix<int> idx(M.m_rows, 1);
  156 + LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr());
  157 + return idx;
  158 + }
  159 + static Matrix<int> getrf(Matrix<double> M) {
  160 + Matrix<int> idx(M.m_rows, 1);
  161 + LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr());
  162 + return idx;
  163 + }
  164 + };
  165 +}
0 166 \ No newline at end of file
... ...