Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 75cf8d04 authored by hhakim's avatar hhakim
Browse files

Add scipy namespace.

parent 3bd15f2d
No related branches found
No related tags found
No related merge requests found
...@@ -36,60 +36,85 @@ ...@@ -36,60 +36,85 @@
#ifndef __FAUST_SCIPY__ #ifndef __FAUST_SCIPY__
#define __FAUST_SCIPY__ #define __FAUST_SCIPY__
namespace scipy
/* {
* Slice columns given as an array of indices (pass 1). /*
* This pass counts idx entries and computes a new indptr. * Slice columns given as an array of indices (pass 1).
* * This pass counts idx entries and computes a new indptr.
* Input Arguments: *
* I n_idx - number of indices to slice * Input Arguments:
* I col_idxs[n_idx] - indices to slice * I n_idx - number of indices to slice
* I n_row - major axis dimension * I col_idxs[n_idx] - indices to slice
* I n_col - minor axis dimension * I n_row - major axis dimension
* I Ap[n_row+1] - indptr * I n_col - minor axis dimension
* I Aj[nnz(A)] - indices * I Ap[n_row+1] - indptr
* * I Aj[nnz(A)] - indices
* Output Arguments: *
* I col_offsets[n_col] - cumsum of index repeats * Output Arguments:
* I Bp[n_row+1] - new indptr * I col_offsets[n_col] - cumsum of index repeats
* * I Bp[n_row+1] - new indptr
*/ *
template<class I, class J> */
void csr_column_index1(const J n_idx, template<class I, class J>
const J col_idxs[], void csr_column_index1(const J n_idx,
const J n_row, const J col_idxs[],
const J n_col, const J n_row,
const I Ap[], const J n_col,
const I Aj[], const I Ap[],
I col_offsets[], const I Aj[],
I Bp[]); I col_offsets[],
I Bp[]);
/* /*
* Slice columns given as an array of indices (pass 2). * Slice columns given as an array of indices (pass 2).
* This pass populates indices/data entries for selected columns. * This pass populates indices/data entries for selected columns.
* *
* Input Arguments: * Input Arguments:
* I col_order[n_idx] - order of col indices * I col_order[n_idx] - order of col indices
* I col_offsets[n_col] - cumsum of col index counts * I col_offsets[n_col] - cumsum of col index counts
* I nnz - nnz(A) * I nnz - nnz(A)
* I Aj[nnz(A)] - column indices * I Aj[nnz(A)] - column indices
* T Ax[nnz(A)] - data * T Ax[nnz(A)] - data
* *
* Output Arguments: * Output Arguments:
* I Bj[nnz(B)] - new column indices * I Bj[nnz(B)] - new column indices
* T Bx[nnz(B)] - new data * T Bx[nnz(B)] - new data
* *
*/ */
template<class I, class J, class T> template<class I, class J, class T>
void csr_column_index2(const J col_order[], void csr_column_index2(const J col_order[],
const I col_offsets[], const I col_offsets[],
const J nnz, const J nnz,
const I Aj[], const I Aj[],
const T Ax[], const T Ax[],
I Bj[], I Bj[],
T Bx[]); T Bx[]);
/*
* Slice rows given as an array of indices.
*
* Input Arguments:
* I n_row_idx - number of row indices
* I rows[n_row_idx] - row indices for indexing
* I Ap[n_row+1] - row pointer
* I Aj[nnz(A)] - column indices
* T Ax[nnz(A)] - data
*
* Output Arguments:
* I Bj - new column indices
* T Bx - new data
*
*/
template<class I, class J, class T>
void csr_row_index(const J n_row_idx,
const J rows[],
const I Ap[],
const I Aj[],
const T Ax[],
I Bj[],
T Bx[]);
}
#include "faust_scipy.hpp" #include "faust_scipy.hpp"
#endif #endif
template<class I, class J> namespace scipy
void csr_column_index1(const J n_idx,
const J col_idxs[],
const J n_row,
const J n_col,
const I Ap[],
const I Aj[],
I col_offsets[],
I Bp[])
{ {
// bincount(col_idxs) template<class I, class J>
for(I jj = 0; jj < n_idx; jj++){ void csr_column_index1(const J n_idx,
const I j = col_idxs[jj]; const J col_idxs[],
col_offsets[j]++; const J n_row,
} const J n_col,
const I Ap[],
const I Aj[],
I col_offsets[],
I Bp[])
{
// bincount(col_idxs)
for(I jj = 0; jj < n_idx; jj++){
const I j = col_idxs[jj];
col_offsets[j]++;
}
// Compute new indptr // Compute new indptr
I new_nnz = 0; I new_nnz = 0;
Bp[0] = 0; Bp[0] = 0;
for(I i = 0; i < n_row; i++){ for(I i = 0; i < n_row; i++){
for(I jj = Ap[i]; jj < Ap[i+1]; jj++){ for(I jj = Ap[i]; jj < Ap[i+1]; jj++){
new_nnz += col_offsets[Aj[jj]]; new_nnz += col_offsets[Aj[jj]];
} }
Bp[i+1] = new_nnz; Bp[i+1] = new_nnz;
} }
// cumsum in-place // cumsum in-place
for(I j = 1; j < n_col; j++){ for(I j = 1; j < n_col; j++){
col_offsets[j] += col_offsets[j - 1]; col_offsets[j] += col_offsets[j - 1];
} }
} }
template<class I, class J, class T>
void csr_column_index2(const J col_order[],
const I col_offsets[],
const J nnz,
const I Aj[],
const T Ax[],
I Bj[],
T Bx[])
{
I n = 0;
for(I jj = 0; jj < nnz; jj++){
const I j = Aj[jj];
const I offset = col_offsets[j];
const I prev_offset = j == 0 ? 0 : col_offsets[j-1];
if (offset != prev_offset) {
const T v = Ax[jj];
for(I k = prev_offset; k < offset; k++){
Bj[n] = col_order[k];
Bx[n] = v;
n++;
}
}
}
}
template<class I, class J, class T>
void csr_column_index2(const J col_order[],
const I col_offsets[],
const J nnz,
const I Aj[],
const T Ax[],
I Bj[],
T Bx[])
{
I n = 0;
for(I jj = 0; jj < nnz; jj++){
const I j = Aj[jj];
const I offset = col_offsets[j];
const I prev_offset = j == 0 ? 0 : col_offsets[j-1];
if (offset != prev_offset) {
const T v = Ax[jj];
for(I k = prev_offset; k < offset; k++){
Bj[n] = col_order[k];
Bx[n] = v;
n++;
}
}
}
}
template<class I, class J, class T>
void csr_row_index(const J n_row_idx,
const J rows[],
const I Ap[],
const I Aj[],
const T Ax[],
I Bj[],
T Bx[])
{
for(J i = 0; i < n_row_idx; i++){
const J row = rows[i];
const I row_start = Ap[row];
const I row_end = Ap[row+1];
Bj = std::copy(Aj + row_start, Aj + row_end, Bj);
Bx = std::copy(Ax + row_start, Ax + row_end, Bx);
}
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment