本文整理汇总了C++中cs_spalloc函数的典型用法代码示例。如果您正苦于以下问题:C++ cs_spalloc函数的具体用法?C++ cs_spalloc怎么用?C++ cs_spalloc使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了cs_spalloc函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: NM_triplet_to_csr
CSparseMatrix* NM_triplet_to_csr(CSparseMatrix* triplet)
{
#ifdef WITH_MKL_SPBLAS
assert(triplet);
CHECK_MKL(load_mkl_function("mkl_dcsrcoo", (void**)&mkl_dcsrcoo_p));
if (triplet->m != triplet->n)
{
fprintf(stderr, "NM_triplet_to_csr :: the matrix has to be square\n");
exit(EXIT_FAILURE);
}
CSparseMatrix* csr = cs_spalloc(NSM_NROW_CSR(triplet), NSM_NCOL_CSR(triplet), triplet->nz, 1, 0);
assert(csr);
csr->nz = -2;
CS_INT n = csr->n;
CS_INT job[6] = {0};
CS_INT info = 0;
job[0] = 2;
(*mkl_dcsrcoo_p)(job, &n, csr->x, csr->i, csr->p, &(triplet->nz), triplet->x, triplet->i, triplet->p, &info);
return csr;
#else
fprintf(stderr, "NM_triplet_to_csr :: MKL not enabled\n");
exit(EXIT_FAILURE);
#endif
}
示例2: tf_calc_dktil
/**
* @brief Creates the penalty matrix D tilde of order k.
* Returns the matrix Dk premultipied by a diagonal
* matrix of weights.
*
* @param n number of observations
* @param k order of the trendfilter
* @param x locations of the responses
* @return pointer to a csparse matrix
* @see tf_calc_dktil
*/
cs * tf_calc_dktil (int n, int k, const double * x)
{
cs * delta_k;
cs * delta_k_cp;
cs * Dk;
cs * Dktil;
int i;
Dk = tf_calc_dk(n, k, x);
/* Deal with k=0 separately */
if(k == 0)
return Dk;
/* Construct diagonal matrix of differences: */
delta_k = cs_spalloc(n-k, n-k, (n-k), 1, 1);
for(i = 0; i < n - k; i++)
{
delta_k->p[i] = i;
delta_k->i[i] = i;
delta_k->x[i] = k / (x[k + i] - x[i]);
}
delta_k->nz = n-k;
delta_k_cp = cs_compress(delta_k);
Dktil = cs_multiply(delta_k_cp, Dk);
cs_spfree(Dk);
cs_spfree(delta_k);
cs_spfree(delta_k_cp);
return Dktil;
}
示例3: NM_csr_to_csc
CSparseMatrix* NM_csr_to_csc(CSparseMatrix* csr)
{
#ifdef WITH_MKL_SPBLAS
assert(csr);
CHECK_MKL(load_mkl_function("mkl_dcsrcsc", (void**)&mkl_dcsrcsc_p));
if (csr->m != csr->n)
{
fprintf(stderr, "NM_csr_to_csc :: the matrix has to be square\n");
exit(EXIT_FAILURE);
}
assert(csr);
CSparseMatrix* csc = cs_spalloc(csr->m, csr->n, csr->nzmax, 1, 0);
assert(csc);
CS_INT n = csr->n;
CS_INT job[6] = {0};
job[5] = 1;
CS_INT info = 0;
(*mkl_dcsrcsc_p)(job, &n, csr->x, csr->i, csr->p, csc->x, csc->i, csc->p, &info);
return csc;
#else
fprintf(stderr, "NM_csr_to_csc :: MKL not enabled\n");
exit(EXIT_FAILURE);
#endif
}
示例4: NM_csr_to_triplet
CSparseMatrix* NM_csr_to_triplet(CSparseMatrix* csr)
{
#ifdef WITH_MKL_SPBLAS
assert(csr);
CHECK_MKL(load_mkl_function("mkl_dcsrcoo", (void**)&mkl_dcsrcoo_p));
if (csr->m != csr->n)
{
fprintf(stderr, "NM_csr_to_triplet :: the matrix has to be square\n");
exit(EXIT_FAILURE);
}
CSparseMatrix* triplet = cs_spalloc(csr->m, csr->n, csr->p[csr->m], 1, 1);
assert(triplet);
CS_INT n = csr->n;
CS_INT job[6] = {0};
job[4] = csr->p[csr->m];
job[5] = 3;
CS_INT info = 0;
(*mkl_dcsrcoo_p)(job, &n, csr->x, csr->i, csr->p, &(csr->p[csr->m]), triplet->x, triplet->i, triplet->p, &info);
triplet->nz = csr->p[csr->m];
return triplet;
#else
fprintf(stderr, "NM_csr_to_triplet :: MKL not enabled\n");
exit(EXIT_FAILURE);
#endif
}
示例5: scalar_plus_diag
/* Calculates A*b + D (identity), for a square matrix
A and scalar b */
cs * scalar_plus_diag (const cs * A, double b, double *D)
{
int i;
int j;
cs * B;
B = cs_spalloc(A->m, A->n, A->nzmax, 1, 0);
for (j = 0; j < A->n; j++)
{
B->p[j] = A->p[j];
for (i = A->p[j] ; i < A->p[j+1] ; i++)
{
if(A->i[i] == j)
{
B->x[i] = b * A->x[i] + D[j];
} else {
B->x[i] = b * A->x[i];
}
B->i[i] = A->i[i];
}
}
B->p[j] = A->p[j];
return B;
}
示例6: return
cs *cs_compress(const cs *T) {
int m, n, nz, p, k, *Cp, *Ci, *w, *Ti, *Tj;
double *Cx, *Tx;
cs *C;
if (!CS_TRIPLET (T))
return (NULL); /* check inputs */
m = T->m;
n = T->n;
Ti = T->i;
Tj = T->p;
Tx = T->x;
nz = T->nz;
C = cs_spalloc(m, n, nz, Tx != NULL, 0); /* allocate result */
w = (int *) cs_calloc(n, sizeof(int)); /* get workspace */
if (!C || !w)
return (cs_done(C, w, NULL, 0)); /* out of memory */
Cp = C->p;
Ci = C->i;
Cx = C->x;
for (k = 0; k < nz; k++)
w[Tj[k]]++; /* column counts */
cs_cumsum(Cp, w, n); /* column pointers */
for (k = 0; k < nz; k++) {
Ci[p = w[Tj[k]]++] = Ti[k]; /* A(i,j) is the pth entry in C */
if (Cx)
Cx[p] = Tx[k];
}
return (cs_done(C, w, NULL, 1)); /* success; free w and return C */
}
示例7: ACADOERROR
returnValue ACADOcsparse::setMatrix(double *A_)
{
int run1;
int order = 0;
if (dim <= 0)
return ACADOERROR(RET_MEMBER_NOT_INITIALISED);
if (nDense <= 0)
return ACADOERROR(RET_MEMBER_NOT_INITIALISED);
cs *C, *D;
C = cs_spalloc(0, 0, 1, 1, 1);
for (run1 = 0; run1 < nDense; run1++)
cs_entry(C, index1[run1], index2[run1], A_[run1]);
D = cs_compress(C);
S = cs_sqr(order, D, 0);
N = cs_lu(D, S, TOL);
cs_spfree(C);
cs_spfree(D);
return SUCCESSFUL_RETURN;
}
示例8: cs_calloc
/* C = A*B */
cs *cs_multiply (const cs *A, const cs *B)
{
CS_INT p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values, *Bi ;
CS_ENTRY *x, *Bx, *Cx ;
cs *C ;
if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ; /* check inputs */
if (A->n != B->m) return (NULL) ;
m = A->m ; anz = A->p [A->n] ;
n = B->n ; Bp = B->p ; Bi = B->i ; Bx = B->x ; bnz = Bp [n] ;
w = cs_calloc (m, sizeof (CS_INT)) ; /* get workspace */
values = (A->x != NULL) && (Bx != NULL) ;
x = values ? cs_malloc (m, sizeof (CS_ENTRY)) : NULL ; /* get workspace */
C = cs_spalloc (m, n, anz + bnz, values, 0) ; /* allocate result */
if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ;
Cp = C->p ;
for (j = 0 ; j < n ; j++)
{
if (nz + m > C->nzmax && !cs_sprealloc (C, 2*(C->nzmax)+m))
{
return (cs_done (C, w, x, 0)) ; /* out of memory */
}
Ci = C->i ; Cx = C->x ; /* C->i and C->x may be reallocated */
Cp [j] = nz ; /* column j of C starts here */
for (p = Bp [j] ; p < Bp [j+1] ; p++)
{
nz = cs_scatter (A, Bi [p], Bx ? Bx [p] : 1, w, x, j+1, C, nz) ;
}
if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ;
}
Cp [n] = nz ; /* finalize the last column of C */
cs_sprealloc (C, 0) ; /* remove extra space from C */
return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */
}
示例9: cs_spalloc
cs *cs_rR(const cs *A, double nu, double nuR, const css *As, const cs *Roldinv, double Roldldet, const cs *pG){
cs *Rnew, *Rnewinv, *Ainv;
double Rnewldet, MH;
int dimG = A->n;
int cnt = 0;
int i, j;
Rnewinv = cs_spalloc (dimG, dimG, dimG*dimG, 1, 0);
for (i = 0 ; i < dimG; i++){
Rnewinv->p[i] = i*dimG;
for (j = 0 ; j < dimG; j++){
Rnewinv->i[cnt] = j;
Rnewinv->x[cnt] = 0.0;
A->x[i*dimG+j] -= pG->x[i*dimG+j];
cnt++;
}
}
Rnewinv->p[dimG] = dimG*dimG;
cs_cov2cor(A);
Ainv = cs_inv(A);
Rnew = cs_rinvwishart(Ainv, nu, As);
cs_cov2cor(Rnew);
Rnewldet = log(cs_invR(Rnew, Rnewinv));
/*****************************************************/
/* From Eq A.4 in Liu and Daniels (2006) */
/* using \pi_{1} = Eq 6 in Barnard (2000) */
/* using \pi_{2} = Eq 3.4 in Liu and Daniels (2006) */
/*****************************************************/
MH = Roldldet-Rnewldet;
for (i = 0 ; i < dimG; i++){
MH += log(Roldinv->x[i*dimG+i]);
MH -= log(Rnewinv->x[i*dimG+i]);
}
MH *= 0.5*nuR;
if(MH<log(runif(0.0,1.0)) || Rnewldet<log(Dtol)){
Rnewldet = cs_invR(Roldinv, Rnew); // save old R
}
for (i = 0 ; i < dimG; i++){
for (j = 0 ; j < dimG; j++){
Rnew->x[i*dimG+j] *= sqrt((pG->x[i*dimG+i])*(pG->x[j*dimG+j]));
}
}
cs_spfree(Rnewinv);
cs_spfree(Ainv);
return (cs_done (Rnew, NULL, NULL, 1)) ; /* success; free workspace, return C */
}
示例10: cs_calloc
/* C = alpha*A + beta*B */
cs *cs_add (const cs *A, const cs *B, double alpha, double beta)
{
int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values ;
double *x, *Bx, *Cx ;
cs *C ;
if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ; /* check inputs */
if (A->m != B->m || A->n != B->n) return (NULL) ;
m = A->m ; anz = A->p [A->n] ;
n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [n] ;
w = cs_calloc (m, sizeof (int)) ; /* get workspace */
values = (A->x != NULL) && (Bx != NULL) ;
x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */
C = cs_spalloc (m, n, anz + bnz, values, 0) ; /* allocate result*/
if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ;
Cp = C->p ; Ci = C->i ; Cx = C->x ;
for (j = 0 ; j < n ; j++)
{
Cp [j] = nz ; /* column j of C starts here */
nz = cs_scatter (A, j, alpha, w, x, j+1, C, nz) ; /* alpha*A(:,j)*/
nz = cs_scatter (B, j, beta, w, x, j+1, C, nz) ; /* beta*B(:,j) */
if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ;
}
Cp [n] = nz ; /* finalize the last column of C */
cs_sprealloc (C, 0) ; /* remove extra space from C */
return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */
}
示例11: ppp_copy
void ppp_copy(cs * to, cs const * from)
{
int j;
if(!CS_CSC(to) || !CS_CSC(from)) REPORT_ERR(PPP_EOP);
if(to->n != from->n || to->m != from->m) REPORT_ERR(PPP_EOP);
cs * zero = cs_spalloc(from->m, from->n, 0, 1, 0);
for(j=0;j<=zero->n;++j)
zero->p[j] = 0;
ppp_add(to, from, zero, 1);
cs_spfree(zero);
}
示例12: econ_createGMatrix
/**
* @brief Creates the admittance matrix.
*
* @return 0 on success.
*/
static int econ_createGMatrix (void)
{
int ret;
int i, j;
double R, Rsum;
cs *M;
StarSystem *sys;
/* Create the matrix. */
M = cs_spalloc( systems_nstack, systems_nstack, 1, 1, 1 );
if (M == NULL)
ERR("Unable to create CSparse Matrix.");
/* Fill the matrix. */
for (i=0; i < systems_nstack; i++) {
sys = &systems_stack[i];
Rsum = 0.;
/* Set some values. */
for (j=0; j < sys->njumps; j++) {
/* Get the resistances. */
R = econ_calcJumpR( sys, &systems_stack[sys->jumps[j]] );
R = 1./R; /* Must be inverted. */
Rsum += R;
/* Matrix is symetrical and non-diagonal is negative. */
ret = cs_entry( M, i, sys->jumps[j], -R );
if (ret != 1)
WARN("Unable to enter CSparse Matrix Cell.");
ret = cs_entry( M, sys->jumps[j], i, -R );
if (ret != 1)
WARN("Unable to enter CSparse Matrix Cell.");
}
/* Set the diagonal. */
Rsum += 1./ECON_SELF_RES; /* We add a resistence for dampening. */
cs_entry( M, i, i, Rsum );
}
/* Compress M matrix and put into G. */
if (econ_G != NULL)
cs_spfree( econ_G );
econ_G = cs_compress( M );
if (econ_G == NULL)
ERR("Unable to create economy G Matrix.");
/* Clean up. */
cs_spfree(M);
return 0;
}
示例13: cs_spalloc
/* load a triplet matrix from a file */
cs *cs_load (FILE *f)
{
int i, j ;
double x ;
cs *T ;
if (!f) return (NULL) ; /* check inputs */
T = cs_spalloc (0, 0, 1, 1, 1) ; /* allocate result */
while (fscanf (f, "%d %d %lg\n", &i, &j, &x) == 3)
{
if (!cs_entry (T, i, j, x)) return (cs_spfree (T)) ;
}
return (T) ;
}
示例14: mexErrMsgTxt
/* cs_updown: sparse Cholesky update/downdate (rank-1 or multiple rank) */
void mexFunction
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
cs Lmatrix, *Lin, Cmatrix, *C, *L, Cvector, *Cvec ;
int ignore, j, k, n, lnz, *parent, sigma = 1, cp [2] ;
char sigma_string [20] ;
if (nargout > 1 || nargin < 3 || nargin > 4)
{
mexErrMsgTxt ("Usage: L = cs_updown(L,C,parent,sigma)") ;
}
Lin = cs_mex_get_sparse (&Lmatrix, 1, 1, pargin [0]) ; /* get input L */
n = Lin->n ;
if (nargin > 3 && mxIsChar (pargin [3]))
{
mxGetString (pargin [3], sigma_string, 8) ;
sigma = (sigma_string [0] == '-') ? (-1) : 1 ;
}
/* make a copy of L (this can take more work than updating L itself) */
lnz = Lin->p [n] ;
L = cs_spalloc (n, n, lnz, 1, 0) ;
for (j = 0 ; j <= n ; j++) L->p [j] = Lin->p [j] ;
for (k = 0 ; k < lnz ; k++) L->i [k] = Lin->i [k] ;
for (k = 0 ; k < lnz ; k++) L->x [k] = Lin->x [k] ;
cs_mex_check (0, n, -1, 0, 1, 1, pargin [1]) ; /* get C */
C = cs_mex_get_sparse (&Cmatrix, 0, 1, pargin [1]) ;
parent = cs_mex_get_int (n, pargin [2], &ignore, 0) ; /* get parent */
/* do the update one column at a time */
Cvec = &Cvector ;
Cvec->m = n ;
Cvec->n = 1 ;
Cvec->p = cp ;
Cvec->nz = -1 ;
cp [0] = 0 ;
for (k = 0 ; k < C->n ; k++)
{
/* extract C(:,k) */
cp [1] = C->p [k+1] - C->p [k] ;
Cvec->nzmax = cp [1] ;
Cvec->i = C->i + C->p [k] ;
Cvec->x = C->x + C->p [k] ;
cs_updown (L, sigma, Cvec, parent) ; /* update/downdate */
}
pargout [0] = cs_mex_put_sparse (&L) ; /* return new L */
}
示例15: formKKT
cs * formKKT(const AMatrix * A, const Settings * s) {
/* ONLY UPPER TRIANGULAR PART IS STUFFED
* forms column compressed KKT matrix
* assumes column compressed form A matrix
*
* forms upper triangular part of [I A'; A -I]
*/
scs_int j, k, kk;
cs * K_cs;
/* I at top left */
const scs_int Anz = A->p[A->n];
const scs_int Knzmax = A->n + A->m + Anz;
cs * K = cs_spalloc(A->m + A->n, A->m + A->n, Knzmax, 1, 1);
#if EXTRAVERBOSE > 0
scs_printf("forming KKT\n");
#endif
if (!K) {
return NULL;
}
kk = 0;
for (k = 0; k < A->n; k++) {
K->i[kk] = k;
K->p[kk] = k;
K->x[kk] = s->rho_x;
kk++;
}
/* A^T at top right : CCS: */
for (j = 0; j < A->n; j++) {
for (k = A->p[j]; k < A->p[j + 1]; k++) {
K->p[kk] = A->i[k] + A->n;
K->i[kk] = j;
K->x[kk] = A->x[k];
kk++;
}
}
/* -I at bottom right */
for (k = 0; k < A->m; k++) {
K->i[kk] = k + A->n;
K->p[kk] = k + A->n;
K->x[kk] = -1;
kk++;
}
/* assert kk == Knzmax */
K->nz = Knzmax;
K_cs = cs_compress(K);
cs_spfree(K);
return (K_cs);
}