Getting Started

The following is a long, running example highlighting most of the features in slap, and is included in the repository at examples/getting_started.c.

  1#include <stdlib.h>
  2#include <math.h>
  3#include "slap/slap.h"
  4
  5sfloat myfun(sfloat x) {
  6  return 2 * x * sin(x);
  7}
  8
  9int main(void) {
 10  puts("Welcome to Getting Started with the slap Library!");
 11
 12  /////////////////////////////////////////////
 13  // Basic Operations
 14  /////////////////////////////////////////////
 15
 16  printf("\n~~~~~~~~~~~~~~~ BASIC OPS ~~~~~~~~~~~~~~\n");
 17
 18  // Create a matrix with stack-allocated memory
 19  sfloat data_A[12] = {1,2,3, 4,5,6, 7,8,9, 10,11,12};
 20  Matrix A = slap_MatrixFromArray(3, 4, data_A);
 21
 22  // Get the sizes of the matrix
 23  int n_rows = slap_NumRows(A);
 24  int n_cols = slap_NumCols(A);
 25  int n_el = slap_NumElements(A);
 26  printf("Size of A = (%d,%d) with %d elements\n", n_rows, n_cols, n_el);
 27
 28  // Get a pointer to an element
 29  sfloat *pa = slap_GetElement(A, 0, 0);
 30  sfloat a = *slap_GetElement(A, 0, 1);
 31  printf("A[0,0] = %0.3g\n", *pa);
 32  printf("A[0,1] = %0.3g\n", a);
 33
 34  // Set an element via pointer
 35  *pa = -20;
 36  printf("A[0,0] = %0.3g (after set)\n", *slap_GetElement(A, 0, 0));
 37
 38  // Set an element directly
 39  slap_SetElement(A, 1, 2, 25.5);
 40  printf("A[1,2] = %0.3g (after set)\n", *slap_GetElement(A, 1, 2));
 41
 42  /////////////////////////////////////////////
 43  // Unary Operations
 44  /////////////////////////////////////////////
 45
 46  printf("\n~~~~~~~~~~~~~~~ UNARY OPS ~~~~~~~~~~~~~~\n");
 47
 48  // Printing
 49  printf("\nA:\n");
 50  slap_PrintMatrix(A);
 51
 52  // Set the matrix to a constant
 53  const sfloat val = 1.2;
 54  slap_SetConst(A, val);
 55  printf("\nA: (set const)\n");
 56  slap_PrintMatrix(A);
 57
 58  // Scale by constant
 59  slap_ScaleByConst(A, -1);
 60  printf("\nA: (Scale by const)\n");
 61  slap_PrintMatrix(A);
 62
 63  // Set Identity
 64  slap_SetIdentity(A, 2.1);
 65  printf("\nA: (Set to identity)\n");
 66  slap_PrintMatrix(A);
 67
 68  // Set Diagonal
 69  sfloat diag[3] = {1,2,3};
 70  slap_SetDiagonal(A, diag, 3);
 71  printf("\nA: (Set Diagonal)\n");
 72  slap_PrintMatrix(A);
 73
 74  // Set first n elements of diagonal
 75  diag[0] = 10;
 76  diag[1] = 11;
 77  slap_SetDiagonal(A, diag, 2);
 78  printf("\nA: (Set Partial Diagonal)\n");
 79  slap_PrintMatrix(A);
 80
 81  // Set Range
 82  slap_SetRange(A, 1, slap_NumElements(A));
 83  printf("\nA: (Set Range)\n");
 84  slap_PrintMatrix(A);
 85
 86  slap_SetRange(A, 0, 1);
 87  printf("\nA: (Set Range 0 to 1)\n");
 88  slap_PrintMatrix(A);
 89
 90
 91  /////////////////////////////////////////////
 92  // Transformations
 93  /////////////////////////////////////////////
 94
 95  printf("\n~~~~~~~~~~~~ TRANSFORMATIONS ~~~~~~~~~~~\n");
 96
 97  // Transpose
 98  Matrix At = slap_Transpose(A);
 99  printf("\nA transpose:\n");
100  slap_PrintMatrix(At);
101  bool tA = slap_IsTransposed(At);
102  printf("At is transposed? %d\n", tA);
103
104  // Flatten
105  Matrix vec_a = slap_Flatten(A);
106  printf("\nA Flatten\n");
107  slap_PrintMatrix(vec_a);
108
109  // Reshape
110  Matrix A2 = slap_Reshape(A, 2, 6);
111  printf("\nA2\n");
112  slap_PrintMatrix(A2);
113
114  // Modify reshape (note that it changes the original)
115  slap_SetElement(A2, 0, 1, 100);
116  printf("\nA (after edit via reshape)\n");
117  slap_PrintMatrix(A);
118
119  // Reshape to smaller
120  // note this is 1st 4 elements, not top left corner
121  Matrix A_resize = slap_Reshape(A, 2, 2);
122  printf("\nA resize\n");
123  slap_PrintMatrix(A_resize);
124
125  /////////////////////////////////////////////
126  // Copying
127  /////////////////////////////////////////////
128
129  printf("\n~~~~~~~~~~~~~~~ COPYING ~~~~~~~~~~~~~~~\n");
130
131  // Matrix with heap-allocated memory
132  sfloat *data_B = (sfloat*)malloc(n_el * sizeof(sfloat));
133  Matrix B = slap_MatrixFromArray(4, 3, data_B);
134
135  // Copy from transposed array
136  slap_Copy(B, At);
137
138  printf("\nCopy from A to B:\n");
139  printf("A:\n");
140  slap_PrintMatrix(A);
141  printf("B:\n");
142  slap_PrintMatrix(B);
143
144  // Copy from array
145  sfloat data_C[4] = {-1,2,-3,4};
146  slap_CopyFromArray(A_resize, data_C);  // note we're copying to a reshaped version of A
147
148  printf("\nA (after array copy):\n");
149  slap_PrintMatrix(A);
150
151  /////////////////////////////////////////////
152  // Sub-Arrays
153  /////////////////////////////////////////////
154
155  printf("\n~~~~~~~~~~~~~~ SUB-ARRAYS ~~~~~~~~~~~~~~\n");
156
157  // Get view of top-left 2x2 corner
158  Matrix A_sub = slap_CreateSubMatrix(A, 0, 0, 2, 2);
159  printf("\nTop-left corner of A\n");
160  slap_PrintMatrix(A_sub);
161  printf("A is Dense?     %d\n", slap_IsDense(A));
162  printf("A sub is Dense? %d\n", slap_IsDense(A_sub));
163
164  // Copy to Sub-matrix
165  data_C[3] = -50;
166  slap_CopyFromArray(A_sub, data_C);
167  printf("\nA (after copying to sub-matrix)\n");
168  slap_PrintMatrix(A);
169
170  // Get middle elements
171  A_sub = slap_CreateSubMatrix(A, 1, 1, 1, 2);
172  printf("\nMiddle of A\n");
173  slap_PrintMatrix(A_sub);
174
175  /////////////////////////////////////////////
176  // Vector Operations
177  /////////////////////////////////////////////
178
179  printf("\n~~~~~~~~~~~~~~ VECTOR OPS ~~~~~~~~~~~~~~\n");
180
181  // Create Matrices (vectors) on the heap
182  // Note the calls to slap_FreeMatrix at the bottom of this function
183  Matrix x = slap_NewMatrix(3, 1);
184  Matrix y = slap_NewMatrix(3, 1);
185  Matrix z = slap_NewMatrixZeros(3, 1);
186
187  // Set some values
188  slap_SetRange(x, 1, 3);
189  slap_SetConst(y, 1.5);
190
191  // Inner product
192  sfloat dot_xy = slap_InnerProduct(x, y);
193  printf("Inner product = %0.3g\n", dot_xy);
194
195  // Cross product
196  slap_CrossProduct(z, x, y);
197  printf("Cross product: ");
198  slap_PrintMatrix(slap_Transpose(z));
199
200  // Norms (only for dense)
201  printf("Two Norm Squared: %0.2g\n", slap_NormTwoSquared(z));
202  printf("Two Norm:         %0.2g\n", slap_NormTwo(z));
203  printf("One Norm:         %0.2g\n", slap_NormOne(z));
204  printf("Inf Norm:         %0.2g\n", slap_NormInf(z));
205
206  // Max/Min
207  printf("Max: %0.2g\n", slap_Max(z));
208  printf("Min: %0.2g\n", slap_Min(z));
209
210  // Argmax/Argmin
211  sfloat z_max, z_min;
212  MatrixIterator it_max = slap_ArgMax(z, &z_max);
213  MatrixIterator it_min = slap_ArgMin(z, &z_min);
214  printf("ArgMax: %0.2g at linear index, %d, Cartesian index (%d,%d)\n", z_max,
215         it_max.k, it_max.i, it_max.j);
216  printf("ArgMin: %0.2g at linear index, %d, Cartesian index (%d,%d)\n", z_max,
217         it_min.k, it_min.i, it_min.j);
218
219  // Outer Product
220  Matrix Z = slap_NewMatrix(3, 3);
221  slap_OuterProduct(Z, x, y);
222  printf("\nOuter Product:\n");
223  slap_PrintMatrix(Z);
224
225  /////////////////////////////////////////////
226  // Linear Algebra
227  /////////////////////////////////////////////
228
229  printf("\n~~~~~~~~~~~~~~ LINEAR ALG ~~~~~~~~~~~~~~\n");
230
231  // Reshape A to be square
232  A = slap_Reshape(A, 3, 3);
233
234  // Create output C matrix
235  Matrix C = slap_NewMatrix(3, 4);
236  slap_SetConst(C, 1);
237
238  // Matrix Multiplication (C = beta * C + alpha * A * B)
239  sfloat alpha = 1.0;
240  sfloat beta = 0.5;
241  slap_MatMulAdd(C, A, slap_Transpose(B), alpha, beta);
242  printf("\nC (MatMulAdd)\n");
243  slap_PrintMatrix(C);
244
245  // Linear Solve w/ Cholesky
246  sfloat data_A2[9];                       // data for PD matrix
247  sfloat data_b[3] = {10.1, -11.2, 12.3};  // rhs vector
248  A2 = slap_MatrixFromArray(3, 3, data_A2);
249  Matrix b = slap_MatrixFromArray(3, 1, data_b);
250
251  slap_MatMulAtB(A2, A, A);    // A2 = A'A. NOTE: Use with caution. Please see docs.
252  slap_AddIdentity(A2, 0.01);  // Ensure A2 > 0
253  slap_Copy(A, A2);      // Save A2 back to A for later
254
255  enum slap_ErrorCode err;
256  err = slap_Cholesky(A2);     // Perform Cholesky decomposition
257  if (err == SLAP_CHOLESKY_FAIL) {
258    printf("Matrix is not Positive-Definite!\n");
259  }
260
261  slap_Copy(x, b);        // copy rhs to x
262  slap_CholeskySolve(A2, x);    // solve for x
263
264  slap_MatMulAB(y, A, x);       // Calculate y = A * x (y should equal b)
265  printf("\nCholesky Solve (these should be equal):\n");
266  printf("b = ");
267  slap_PrintMatrix(slap_Transpose(b));
268  printf("y = ");
269  slap_PrintMatrix(slap_Transpose(y));
270
271  /////////////////////////////////////////////
272  // Advanced Ops
273  /////////////////////////////////////////////
274
275  printf("\n~~~~~~~~~~~~~ ADVANCED OPS ~~~~~~~~~~~~~\n");
276
277  // Efficient Iteration (including Sub-Arrays)
278  Matrix B_sub = slap_CreateSubMatrix(B, 1, 1, 3, 2);
279  printf("\nB\n");
280  slap_PrintMatrix(B);
281
282  printf("\nB_sub\n");
283  slap_PrintMatrix(B_sub);
284
285  printf("\nIteration over Sub-Array:\n");
286  for (MatrixIterator it = slap_Iterator(B_sub); !slap_IsFinished(&it); slap_Step(&it)) {
287    int mem_index = it.index;
288    int lin_index = it.k;
289    int row_index = it.i;
290    int col_index = it.j;
291    sfloat B_val = B_sub.data[mem_index];
292    printf("B_sub[%d,%d] = % 4.2f at linear index %d and memory index %d\n", row_index,
293           col_index, B_val, lin_index, mem_index);
294  }
295
296  // Function Mapping
297  printf("\nB (before map)\n");
298  slap_PrintMatrix(B);
299
300  slap_Map(B_sub, myfun);
301
302  printf("\nB (after map)\n");
303  slap_PrintMatrix(B);
304
305  /////////////////////////////////////////////
306  // Memory Cleanup (DON'T FORGET THIS)
307  /////////////////////////////////////////////
308
309  free(data_B);
310  slap_FreeMatrix(&x);
311  slap_FreeMatrix(&y);
312  slap_FreeMatrix(&z);
313  slap_FreeMatrix(&Z);
314  slap_FreeMatrix(&C);
315  return EXIT_SUCCESS;
316}