Re: Help with [] for defined types

[David Bateman]

Hi Andy,


According to Andy Adler <address@hidden> (on 08/02/04):
> I've done this, and it seems to be being called twice, the first time
> with dv=[0,0] and the second time with the right output size, but
> 
> echo 's=sparse(diag([1,3],1)); p=[s,s]' | octave -qfH
> DEBUG:sparse( void)
> installing sparse type at type-id = 40
> install sm_sm
> DEBUG:sparse - matrix_to_sparse
> DEBUG:sparse( SuperMatrix A)
> DEBUG:sparse_resize                  <--- calls sparse([1;2],[2;3],[1;3],0,0)
> DEBUG:sparse - assemble_sparse
> error: sparse row index out of range  <-- error because rows,cols= 0,0
> DEBUG:sparse( SuperMatrix A)
> DEBUG:sparse_resize                  <--- calls sparse([],[],[],3,6)
> DEBUG:sparse( SuperMatrix A)
> DEBUG:sparse destructor
> DEBUG:sparse destructor
> error: evaluating assignment expression near line 1, column 28
> 
> I don't understand why octave_sparse::resize is being called twice here
> Am I missing something?

The explanation is in the comment in tree_matrix::rvalue, which states

      // The line below might seem crazy, since we take a copy
      // of the first argument, resize it to be empty and then resize
      // it to be full. This is done since it means that there is no
      // recopying of data, as would happen if we used a single resize.
      // It should be noted that resize operation is also significantly 
      // slower than the do_cat_op function, so it makes sense to have an
      // empty matrix and copy all data.
      //
      // We might also start with a empty octave_value using
      //    ctmp = octave_value_typeinfo::lookup_type
      //          (tmp.begin() -> begin() -> type_name());
      // and then directly resize. However, for some types there might be
      // some additional setup needed, and so this should be avoided.

The error comes about because, resize is supposed to be destructive. So
you should modifiy your ASSEMBLE_SPARSE macro to be something like

#define ASSEMBLE_SPARSE2( TYPX, RESIZE) \
   int  nnz= MAX( ridxA.length(), cidxA.length() ); \
   TYPX* coefX= (TYPX*)oct_sparse_malloc( nnz  * sizeof(TYPX)); \
   int * ridxX= (int *)oct_sparse_malloc( nnz  * sizeof(int) ); \
   int * cidxX= (int *)oct_sparse_malloc( (n+1)* sizeof(int)); cidxX[0]= 0; \
 \
   bool ri_scalar = (ridxA.length() == 1); \
   bool ci_scalar = (cidxA.length() == 1); \
   bool cf_scalar = (coefA.length() == 1); \
 \
   sort_idxl* sidx = (sort_idxl *)oct_sparse_malloc( nnz* sizeof(sort_idxl) );\
   int actual_nnz=0; \
   OCTAVE_QUIT; \
   for (int i=0; i<nnz; i++) { \
      double rowidx=  ( ri_scalar ? ridxA(0) : ridxA(i) ) - 1; \
      double colidx=  ( ci_scalar ? cidxA(0) : cidxA(i) ) - 1;  \
      if (rowidx >= m || rowidx < 0) \
         if (! RESIZE ) \
            error("sparse row index out of range"); \
      else if (colidx >= n || colidx < 0) \
         if (! RESIZE ) \
            error("sparse column index out of range"); \
      if (error_state) { actual_nnz=0; break; } \
      if ( coefA(cf_scalar?0:i) != 0. ) { \
         sidx[actual_nnz].val = (long) ( rowidx + m*colidx ); \
         sidx[actual_nnz].idx = i; \
         actual_nnz++; \
      } \
   } \
<snip the rest of the old ASSEMBLE_SPARSE>

#define ASSEMBLE_SPARSE( TYPX  )  ASSEMBLE_SPARSE2 (TYPX, false)

You can then call ASSEMBLE_SPARSE2 in the resize operations

> I've implemented these, but they don't seem to be being called.
> Definitely, octave_sparse concat (const octave_sparse& ra,
>                       const octave_sparse& rb, const Array<int>& ra_idx);
> is not called.
> 
> Also, I'm not sure I understand what ra_idx is and how I'm supposed to
> interpret it.

As the failing resize operations in tree_matrix::rvalue is followed be

      if (error_state)
        goto done;

This is not surprising. Fix up resize so that it is destructive and you'll
probably get to the code that calls concat. Then another set of bugs will
probably appear :-)

Cheers
David

-- 
David Bateman                                address@hidden
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary