sprow.c

Meschach 可以解稠密或稀疏线性方程组、计算特征值和特征向量和解最小平方问题

   len1 = r1->len;	len2 = r2->len;	len_out = r_out->maxlen;
   idx1 = idx2 = idx_out = 0;
   elt1 = r1->elt;	elt2 = r2->elt;	elt_out = r_out->elt;
   
   while ( idx1 < len1 || idx2 < len2 )
   {
      while ( idx_out >= len_out )
      {   /* r_out is too small */
	 r_out->len = idx_out;
	 r_out = sprow_xpd(r_out,0,type);
	 len_out = r_out->maxlen;
	 elt_out = &(r_out->elt[idx_out]);
      }
      if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
      {
	 elt_out->col = elt1->col;
	 elt_out->val = elt1->val;
	 if ( elt1->col == elt2->col && idx2 < len2 )
	 {	elt2++;		idx2++;	}
	 elt1++;	idx1++;
      }
      else
      {
	 elt_out->col = elt2->col;
	 elt_out->val = 0.0;
	 elt2++;	idx2++;
      }
      elt_out++;	idx_out++;
   }
   r_out->len = idx_out;
   
   return r_out;
}

/* sprow_mltadd -- sets r_out <- r1 + alpha.r2
   -- cannot be in situ
   -- only for columns j0, j0+1, ...
   -- type must be TYPE_SPMAT or TYPE_SPROW depending on
      whether r_out is a row of a SPMAT structure
      or a SPROW variable
   -- returns r_out */
#ifndef ANSI_C
SPROW	*sprow_mltadd(r1,r2,alpha,j0,r_out,type)
SPROW	*r1, *r2, *r_out;
double	alpha;
int	j0, type;
#else
SPROW	*sprow_mltadd(const SPROW *r1,const SPROW *r2, double alpha,
		      int j0, SPROW *r_out, int type)
#endif
{
   int	idx1, idx2, idx_out, len1, len2, len_out;
   row_elt	*elt1, *elt2, *elt_out;
   
   if ( ! r1 || ! r2 )
     error(E_NULL,"sprow_mltadd");
   if ( r1 == r_out || r2 == r_out )
     error(E_INSITU,"sprow_mltadd");
   if ( j0 < 0 )
     error(E_BOUNDS,"sprow_mltadd");
   if ( ! r_out )
     r_out = sprow_get(MINROWLEN);
   
   /* Initialise */
   len1 = r1->len;	len2 = r2->len;	len_out = r_out->maxlen;
   /* idx1 = idx2 = idx_out = 0; */
   idx1    = sprow_idx(r1,j0);
   idx2    = sprow_idx(r2,j0);
   idx_out = sprow_idx(r_out,j0);
   idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
   idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
   idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
   elt1    = &(r1->elt[idx1]);
   elt2    = &(r2->elt[idx2]);
   elt_out = &(r_out->elt[idx_out]);
   
   while ( idx1 < len1 || idx2 < len2 )
   {
      if ( idx_out >= len_out )
      {   /* r_out is too small */
	 r_out->len = idx_out;
	 r_out = sprow_xpd(r_out,0,type);
	 len_out = r_out->maxlen;
	 elt_out = &(r_out->elt[idx_out]);
      }
      if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
      {
	 elt_out->col = elt1->col;
	 elt_out->val = elt1->val;
	 if ( idx2 < len2 && elt1->col == elt2->col )
	 {
	    elt_out->val += alpha*elt2->val;
	    elt2++;		idx2++;
	 }
	 elt1++;	idx1++;
      }
      else
      {
	 elt_out->col = elt2->col;
	 elt_out->val = alpha*elt2->val;
	 elt2++;	idx2++;
      }
      elt_out++;	idx_out++;
   }
   r_out->len = idx_out;
   
   return r_out;
}

/* sprow_add -- sets r_out <- r1 + r2
   -- cannot be in situ
   -- only for columns j0, j0+1, ...
   -- type must be TYPE_SPMAT or TYPE_SPROW depending on
      whether r_out is a row of a SPMAT structure
      or a SPROW variable
   -- returns r_out */
#ifndef ANSI_C
SPROW	*sprow_add(r1,r2,j0,r_out,type)
SPROW	*r1, *r2, *r_out;
int	j0, type;
#else
SPROW	*sprow_add(const SPROW *r1,const SPROW *r2, 
		   int j0, SPROW *r_out, int type)
#endif
{
   int	idx1, idx2, idx_out, len1, len2, len_out;
   row_elt	*elt1, *elt2, *elt_out;
   
   if ( ! r1 || ! r2 )
     error(E_NULL,"sprow_add");
   if ( r1 == r_out || r2 == r_out )
     error(E_INSITU,"sprow_add");
   if ( j0 < 0 )
     error(E_BOUNDS,"sprow_add");
   if ( ! r_out )
     r_out = sprow_get(MINROWLEN);
   
   /* Initialise */
   len1 = r1->len;	len2 = r2->len;	len_out = r_out->maxlen;
   /* idx1 = idx2 = idx_out = 0; */
   idx1    = sprow_idx(r1,j0);
   idx2    = sprow_idx(r2,j0);
   idx_out = sprow_idx(r_out,j0);
   idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
   idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
   idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
   elt1    = &(r1->elt[idx1]);
   elt2    = &(r2->elt[idx2]);
   elt_out = &(r_out->elt[idx_out]);
   
   while ( idx1 < len1 || idx2 < len2 )
   {
      if ( idx_out >= len_out )
      {   /* r_out is too small */
	 r_out->len = idx_out;
	 r_out = sprow_xpd(r_out,0,type);
	 len_out = r_out->maxlen;
	 elt_out = &(r_out->elt[idx_out]);
      }
      if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
      {
	 elt_out->col = elt1->col;
	 elt_out->val = elt1->val;
	 if ( idx2 < len2 && elt1->col == elt2->col )
	 {
	    elt_out->val += elt2->val;
	    elt2++;		idx2++;
	 }
	 elt1++;	idx1++;
      }
      else
      {
	 elt_out->col = elt2->col;
	 elt_out->val = elt2->val;
	 elt2++;	idx2++;
      }
      elt_out++;	idx_out++;
   }
   r_out->len = idx_out;
   
   return r_out;
}

/* sprow_sub -- sets r_out <- r1 - r2
   -- cannot be in situ
   -- only for columns j0, j0+1, ...
   -- type must be TYPE_SPMAT or TYPE_SPROW depending on
      whether r_out is a row of a SPMAT structure
      or a SPROW variable
   -- returns r_out */
#ifndef ANSI_C
SPROW	*sprow_sub(r1,r2,j0,r_out,type)
SPROW	*r1, *r2, *r_out;
int	j0, type;
#else
SPROW	*sprow_sub(const SPROW *r1, const SPROW *r2,
		   int j0, SPROW *r_out, int type)
#endif
{
   int	idx1, idx2, idx_out, len1, len2, len_out;
   row_elt	*elt1, *elt2, *elt_out;
   
   if ( ! r1 || ! r2 )
     error(E_NULL,"sprow_sub");
   if ( r1 == r_out || r2 == r_out )
     error(E_INSITU,"sprow_sub");
   if ( j0 < 0 )
     error(E_BOUNDS,"sprow_sub");
   if ( ! r_out )
     r_out = sprow_get(MINROWLEN);
   
   /* Initialise */
   len1 = r1->len;	len2 = r2->len;	len_out = r_out->maxlen;
   /* idx1 = idx2 = idx_out = 0; */
   idx1    = sprow_idx(r1,j0);
   idx2    = sprow_idx(r2,j0);
   idx_out = sprow_idx(r_out,j0);
   idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
   idx2    = (idx2 < 0) ? -(idx2+2) : idx2;
   idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
   elt1    = &(r1->elt[idx1]);
   elt2    = &(r2->elt[idx2]);
   elt_out = &(r_out->elt[idx_out]);
   
   while ( idx1 < len1 || idx2 < len2 )
   {
      if ( idx_out >= len_out )
      {   /* r_out is too small */
	 r_out->len = idx_out;
	 r_out = sprow_xpd(r_out,0,type);
	 len_out = r_out->maxlen;
	 elt_out = &(r_out->elt[idx_out]);
      }
      if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
      {
	 elt_out->col = elt1->col;
	 elt_out->val = elt1->val;
	 if ( idx2 < len2 && elt1->col == elt2->col )
	 {
	    elt_out->val -= elt2->val;
	    elt2++;		idx2++;
	 }
	 elt1++;	idx1++;
      }
      else
      {
	 elt_out->col = elt2->col;
	 elt_out->val = -elt2->val;
	 elt2++;	idx2++;
      }
      elt_out++;	idx_out++;
   }
   r_out->len = idx_out;
   
   return r_out;
}


/* sprow_smlt -- sets r_out <- alpha*r1 
   -- can be in situ
   -- only for columns j0, j0+1, ...
   -- returns r_out */
#ifndef ANSI_C
SPROW	*sprow_smlt(r1,alpha,j0,r_out,type)
SPROW	*r1, *r_out;
double	alpha;
int	j0, type;
#else
SPROW	*sprow_smlt(const SPROW *r1, double alpha, int j0, SPROW *r_out, int type)
#endif
{
   int	idx1, idx_out, len1;
   row_elt	*elt1, *elt_out;
   
   if ( ! r1 )
     error(E_NULL,"sprow_smlt");
   if ( j0 < 0 )
     error(E_BOUNDS,"sprow_smlt");
   if ( ! r_out )
     r_out = sprow_get(MINROWLEN);
   
   /* Initialise */
   len1 = r1->len;
   idx1    = sprow_idx(r1,j0);
   idx_out = sprow_idx(r_out,j0);
   idx1    = (idx1 < 0) ? -(idx1+2) : idx1;
   idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
   elt1    = &(r1->elt[idx1]);

   r_out = sprow_resize(r_out,idx_out+len1-idx1,type);  
   elt_out = &(r_out->elt[idx_out]);

   for ( ; idx1 < len1; elt1++,elt_out++,idx1++,idx_out++ )
   {
      elt_out->col = elt1->col;
      elt_out->val = alpha*elt1->val;
   }

   r_out->len = idx_out;

   return r_out;
}

#ifndef MEX
/* sprow_foutput -- print a representation of r on stream fp */
#ifndef ANSI_C
void	sprow_foutput(fp,r)
FILE	*fp;
SPROW	*r;
#else
void	sprow_foutput(FILE *fp, const SPROW *r)
#endif
{
   int	i, len;
   row_elt	*e;
   
   if ( ! r )
   {
      fprintf(fp,"SparseRow: **** NULL ****\n");
      return;
   }
   len = r->len;
   fprintf(fp,"SparseRow: length: %d\n",len);
   for ( i = 0, e = r->elt; i < len; i++, e++ )
     fprintf(fp,"Column %d: %g, next row: %d, next index %d\n",
	     e->col, e->val, e->nxt_row, e->nxt_idx);
}
#endif


/* sprow_set_val -- sets the j-th column entry of the sparse row r
   -- Note: destroys the usual column & row access paths */
#ifndef ANSI_C
double  sprow_set_val(r,j,val)
SPROW   *r;
int     j;
double  val;
#else
double  sprow_set_val(SPROW *r, int j, double val)
#endif
{
   int  idx, idx2, new_len;
   
   if ( ! r )
     error(E_NULL,"sprow_set_val");
   
   idx = sprow_idx(r,j);
   if ( idx >= 0 )
   {    r->elt[idx].val = val;  return val;     }
   /* else */ if ( idx < -1 )
   {
      /* shift & insert new value */
      idx = -(idx+2);   /* this is the intended insertion index */
      if ( r->len >= r->maxlen )
      {
         r->len = r->maxlen;
         new_len = max(2*r->maxlen+1,5);
         if (mem_info_is_on()) {
            mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),
                        new_len*sizeof(row_elt)); 
         }
         
         r->elt = RENEW(r->elt,new_len,row_elt);
         if ( ! r->elt )        /* can't allocate */
           error(E_MEM,"sprow_set_val");
         r->maxlen = 2*r->maxlen+1;
      }
      for ( idx2 = r->len-1; idx2 >= idx; idx2-- )
        MEM_COPY((char *)(&(r->elt[idx2])),
                 (char *)(&(r->elt[idx2+1])),sizeof(row_elt));
      /************************************************************
        if ( idx < r->len )
        MEM_COPY((char *)(&(r->elt[idx])),(char *)(&(r->elt[idx+1])),
        (r->len-idx)*sizeof(row_elt));
        ************************************************************/
      r->len++;
      r->elt[idx].col = j;
      r->elt[idx].nxt_row = -1;
      r->elt[idx].nxt_idx = -1;
      return r->elt[idx].val = val;
   }
   /* else -- idx == -1, error in index/matrix! */
   return 0.0;
}