📄 sparse.c
字号:
if (!in_situ) sp_zero(C);
}
if (tmp == (SPROW *)NULL && in_situ) {
tmp = sprow_get(MINROWLEN);
MEM_STAT_REG(tmp,TYPE_SPROW);
}
if (in_situ)
for (i=0; i < A->m; i++) {
rc = &(C->row[i]);
sprow_sub(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
sprow_resize(rc,tmp->len,TYPE_SPMAT);
MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
rc->len = tmp->len;
}
else
for (i=0; i < A->m; i++) {
sprow_sub(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
}
C->flag_col = C->flag_diag = FALSE;
#ifdef THREADSAFE
sprow_free(tmp);
#endif
return C;
}
/* C = A+alpha*B, cannot be in situ */
#ifndef ANSI_C
SPMAT *sp_mltadd(A,B,alpha,C)
SPMAT *A, *B, *C;
double alpha;
#else
SPMAT *sp_mltadd(const SPMAT *A, const SPMAT *B, double alpha, SPMAT *C)
#endif
{
int i, in_situ;
SPROW *rc;
STATIC SPROW *tmp = NULL;
if ( ! A || ! B )
error(E_NULL,"sp_mltadd");
if ( A->m != B->m || A->n != B->n )
error(E_SIZES,"sp_mltadd");
if (C == A || C == B)
in_situ = TRUE;
else in_situ = FALSE;
if ( ! C )
C = sp_get(A->m,A->n,5);
else {
if ( C->m != A->m || C->n != A->n )
error(E_SIZES,"sp_mltadd");
if (!in_situ) sp_zero(C);
}
if (tmp == (SPROW *)NULL && in_situ) {
tmp = sprow_get(MINROWLEN);
MEM_STAT_REG(tmp,TYPE_SPROW);
}
if (in_situ)
for (i=0; i < A->m; i++) {
rc = &(C->row[i]);
sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,tmp,TYPE_SPROW);
sprow_resize(rc,tmp->len,TYPE_SPMAT);
MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
rc->len = tmp->len;
}
else
for (i=0; i < A->m; i++) {
sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,
&(C->row[i]),TYPE_SPMAT);
}
C->flag_col = C->flag_diag = FALSE;
#ifdef THREADSAFE
sprow_free(tmp);
#endif
return C;
}
/* B = alpha*A, can be in situ */
#ifndef ANSI_C
SPMAT *sp_smlt(A,alpha,B)
SPMAT *A, *B;
double alpha;
#else
SPMAT *sp_smlt(const SPMAT *A, double alpha, SPMAT *B)
#endif
{
int i;
if ( ! A )
error(E_NULL,"sp_smlt");
if ( ! B )
B = sp_get(A->m,A->n,5);
else
if ( A->m != B->m || A->n != B->n )
error(E_SIZES,"sp_smlt");
for (i=0; i < A->m; i++) {
sprow_smlt(&(A->row[i]),alpha,0,&(B->row[i]),TYPE_SPMAT);
}
return B;
}
/* sp_zero -- zero all the (represented) elements of a sparse matrix */
#ifndef ANSI_C
SPMAT *sp_zero(A)
SPMAT *A;
#else
SPMAT *sp_zero(SPMAT *A)
#endif
{
int i, idx, len;
row_elt *elt;
if ( ! A )
error(E_NULL,"sp_zero");
for ( i = 0; i < A->m; i++ )
{
elt = A->row[i].elt;
len = A->row[i].len;
for ( idx = 0; idx < len; idx++ )
(*elt++).val = 0.0;
}
return A;
}
/* sp_copy2 -- copy sparse matrix (type 2)
-- keeps structure of the OUT matrix */
#ifndef ANSI_C
SPMAT *sp_copy2(A,OUT)
SPMAT *A, *OUT;
#else
SPMAT *sp_copy2(const SPMAT *A, SPMAT *OUT)
#endif
{
int i /* , idx, len1, len2 */;
SPROW *r1, *r2;
STATIC SPROW *scratch = (SPROW *)NULL;
/* row_elt *e1, *e2; */
if ( ! A )
error(E_NULL,"sp_copy2");
if ( ! OUT )
OUT = sp_get(A->m,A->n,10);
if ( ! scratch ) {
scratch = sprow_xpd(scratch,MINROWLEN,TYPE_SPROW);
MEM_STAT_REG(scratch,TYPE_SPROW);
}
if ( OUT->m < A->m )
{
if (mem_info_is_on()) {
mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
A->m*sizeof(SPROW));
}
OUT->row = RENEW(OUT->row,A->m,SPROW);
if ( ! OUT->row )
error(E_MEM,"sp_copy2");
for ( i = OUT->m; i < A->m; i++ )
{
OUT->row[i].elt = NEW_A(MINROWLEN,row_elt);
if ( ! OUT->row[i].elt )
error(E_MEM,"sp_copy2");
else if (mem_info_is_on()) {
mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
}
OUT->row[i].maxlen = MINROWLEN;
OUT->row[i].len = 0;
}
OUT->m = A->m;
}
OUT->flag_col = OUT->flag_diag = FALSE;
/* sp_zero(OUT); */
for ( i = 0; i < A->m; i++ )
{
r1 = &(A->row[i]); r2 = &(OUT->row[i]);
sprow_copy(r1,r2,scratch,TYPE_SPROW);
if ( r2->maxlen < scratch->len )
sprow_xpd(r2,scratch->len,TYPE_SPMAT);
MEM_COPY((char *)(scratch->elt),(char *)(r2->elt),
scratch->len*sizeof(row_elt));
r2->len = scratch->len;
/*******************************************************
e1 = r1->elt; e2 = r2->elt;
len1 = r1->len; len2 = r2->len;
for ( idx = 0; idx < len2; idx++, e2++ )
e2->val = 0.0;
for ( idx = 0; idx < len1; idx++, e1++ )
sprow_set_val(r2,e1->col,e1->val);
*******************************************************/
}
sp_col_access(OUT);
#ifdef THREADSAFE
sprow_free(scratch);
#endif
return OUT;
}
/* sp_resize -- resize a sparse matrix
-- don't destroying any contents if possible
-- returns resized matrix */
#ifndef ANSI_C
SPMAT *sp_resize(A,m,n)
SPMAT *A;
int m, n;
#else
SPMAT *sp_resize(SPMAT *A, int m, int n)
#endif
{
int i, len;
SPROW *r;
if (m < 0 || n < 0)
error(E_NEG,"sp_resize");
if ( ! A )
return sp_get(m,n,10);
if (m == A->m && n == A->n)
return A;
if ( m <= A->max_m )
{
for ( i = A->m; i < m; i++ )
A->row[i].len = 0;
A->m = m;
}
else
{
if (mem_info_is_on()) {
mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
m*sizeof(SPROW));
}
A->row = RENEW(A->row,(unsigned)m,SPROW);
if ( ! A->row )
error(E_MEM,"sp_resize");
for ( i = A->m; i < m; i++ )
{
if ( ! (A->row[i].elt = NEW_A(MINROWLEN,row_elt)) )
error(E_MEM,"sp_resize");
else if (mem_info_is_on()) {
mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
}
A->row[i].len = 0; A->row[i].maxlen = MINROWLEN;
}
A->m = A->max_m = m;
}
/* update number of rows */
A->n = n;
/* do we need to increase the size of start_idx[] and start_row[] ? */
if ( n > A->max_n )
{ /* only have to update the start_idx & start_row arrays */
if (mem_info_is_on())
{
mem_bytes(TYPE_SPMAT,2*A->max_n*sizeof(int),
2*n*sizeof(int));
}
A->start_row = RENEW(A->start_row,(unsigned)n,int);
A->start_idx = RENEW(A->start_idx,(unsigned)n,int);
if ( ! A->start_row || ! A->start_idx )
error(E_MEM,"sp_resize");
A->max_n = n; /* ...and update max_n */
return A;
}
if ( n <= A->n )
/* make sure that all rows are truncated just before column n */
for ( i = 0; i < A->m; i++ )
{
r = &(A->row[i]);
len = sprow_idx(r,n);
if ( len < 0 )
len = -(len+2);
if ( len < 0 )
error(E_MEM,"sp_resize");
r->len = len;
}
return A;
}
/* sp_compact -- removes zeros and near-zeros from a sparse matrix */
#ifndef ANSI_C
SPMAT *sp_compact(A,tol)
SPMAT *A;
double tol;
#else
SPMAT *sp_compact(SPMAT *A, double tol)
#endif
{
int i, idx1, idx2;
SPROW *r;
row_elt *elt1, *elt2;
if ( ! A )
error(E_NULL,"sp_compact");
if ( tol < 0.0 )
error(E_RANGE,"sp_compact");
A->flag_col = A->flag_diag = FALSE;
for ( i = 0; i < A->m; i++ )
{
r = &(A->row[i]);
elt1 = elt2 = r->elt;
idx1 = idx2 = 0;
while ( idx1 < r->len )
{
/* printf("# sp_compact: idx1 = %d, idx2 = %d\n",idx1,idx2); */
if ( fabs(elt1->val) <= tol )
{ idx1++; elt1++; continue; }
if ( elt1 != elt2 )
MEM_COPY(elt1,elt2,sizeof(row_elt));
idx1++; elt1++;
idx2++; elt2++;
}
r->len = idx2;
}
return A;
}
/* sp_mlt (C) Copyright David Stewart and Fabrizio Novalis <novalis@mars.elet.polimi.it> */
/* sp_mlt -- computes out = A*B and returns out */
SPMAT *sp_mlt(const SPMAT *A, const SPMAT *B, SPMAT *out)
{
int i, j, k, idx, cp;
SPROW *rA, *rB, *rout, *rtemp;
double valA;
if ( ! A || ! B )
error(E_NULL,"sp_mlt");
if ( A->n != B->m )
error(E_SIZES,"sp_mlt");
out = sp_resize(out,A->m,B->n);
sp_zero(out);
rtemp = sprow_get(B->n);
for ( i = 0; i < A->m; i++ ) /* per ogni riga */
{
rtemp = sprow_resize(rtemp,0,TYPE_SPROW);
rA = &(A->row[i]);
rout = &(out->row[i]);
for ( idx = 0; idx < rA->len; idx++ ) /* per ogni elemento != 0
della riga corrente */
{
j = rA->elt[idx].col;
valA = rA->elt[idx].val;
rB = &(B->row[j]);
sprow_mltadd(rtemp,rB,valA,0,rout,TYPE_SPMAT);
for ( cp = 0; cp < rout->len; cp++ )
{
rtemp->elt[cp].col = rout->elt[cp].col;
rtemp->elt[cp].val = rout->elt[cp].val;
}
rtemp->len=rout->len;
}
}
return out;
}
/* varying number of arguments */
#ifdef ANSI_C
/* To allocate memory to many arguments.
The function should be called:
sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
where
int m,n,deg;
SPMAT *x, *y, *z,...;
The last argument should be NULL !
m x n is the dimension of matrices x,y,z,...
returned value is equal to the number of allocated variables
*/
int sp_get_vars(int m,int n,int deg,...)
{
va_list ap;
int i=0;
SPMAT **par;
va_start(ap, deg);
while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
*par = sp_get(m,n,deg);
i++;
}
va_end(ap);
return i;
}
/* To resize memory for many arguments.
The function should be called:
sp_resize_vars(m,n,&x,&y,&z,...,NULL);
where
int m,n;
SPMAT *x, *y, *z,...;
The last argument should be NULL !
m X n is the resized dimension of matrices x,y,z,...
returned value is equal to the number of allocated variables.
If one of x,y,z,.. arguments is NULL then memory is allocated to this
argument.
*/
int sp_resize_vars(int m,int n,...)
{
va_list ap;
int i=0;
SPMAT **par;
va_start(ap, n);
while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
*par = sp_resize(*par,m,n);
i++;
}
va_end(ap);
return i;
}
/* To deallocate memory for many arguments.
The function should be called:
sp_free_vars(&x,&y,&z,...,NULL);
where
SPMAT *x, *y, *z,...;
The last argument should be NULL !
There must be at least one not NULL argument.
returned value is equal to the number of allocated variables.
Returned value of x,y,z,.. is VNULL.
*/
int sp_free_vars(SPMAT **va,...)
{
va_list ap;
int i=1;
SPMAT **par;
sp_free(*va);
*va = (SPMAT *) NULL;
va_start(ap, va);
while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
sp_free(*par);
*par = (SPMAT *)NULL;
i++;
}
va_end(ap);
return i;
}
#elif VARARGS
/* To allocate memory to many arguments.
The function should be called:
sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
where
int m,n,deg;
SPMAT *x, *y, *z,...;
The last argument should be NULL !
m x n is the dimension of matrices x,y,z,...
returned value is equal to the number of allocated variables
*/
int sp_get_vars(va_alist) va_dcl
{
va_list ap;
int i=0, m, n, deg;
SPMAT **par;
va_start(ap);
m = va_arg(ap,int);
n = va_arg(ap,int);
deg = va_arg(ap,int);
while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
*par = sp_get(m,n,deg);
i++;
}
va_end(ap);
return i;
}
/* To resize memory for many arguments.
The function should be called:
sp_resize_vars(m,n,&x,&y,&z,...,NULL);
where
int m,n;
SPMAT *x, *y, *z,...;
The last argument should be NULL !
m X n is the resized dimension of matrices x,y,z,...
returned value is equal to the number of allocated variables.
If one of x,y,z,.. arguments is NULL then memory is allocated to this
argument.
*/
int sp_resize_vars(va_alist) va_dcl
{
va_list ap;
int i=0, m, n;
SPMAT **par;
va_start(ap);
m = va_arg(ap,int);
n = va_arg(ap,int);
while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
*par = sp_resize(*par,m,n);
i++;
}
va_end(ap);
return i;
}
/* To deallocate memory for many arguments.
The function should be called:
sp_free_vars(&x,&y,&z,...,NULL);
where
SPMAT *x, *y, *z,...;
The last argument should be NULL !
There must be at least one not NULL argument.
returned value is equal to the number of allocated variables.
Returned value of x,y,z,.. is VNULL.
*/
int sp_free_vars(va_alist) va_dcl
{
va_list ap;
int i=0;
SPMAT **par;
va_start(ap);
while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
sp_free(*par);
*par = (SPMAT *)NULL;
i++;
}
va_end(ap);
return i;
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -