/*
 * call-seq:
 *   diagonal(array)
 *   diagonal(vector)
 *   diagonal(n, scalar)
 *   diagonal(n) { |i| ... }
 *
 * In the first form, create a diagonal matrix with the elements of
 * +array+ along the diagonal.  The second form is like the first but
 * using elements of the column vector +vector+.  In the third form,
 * create an +n+ x +n+ diagonal matrix with diagonal elements of value
 * +scalar+.  In the fourth form, each element <tt>m[i,i]</tt> gets
 * the block result.
 *
 */
VALUE rb_dmatrix_s_diagonal(int argc, VALUE* argv, VALUE klass)
{
    DMatrix* a ;
    VALUE ra ;
    VALUE rb ;
    VALUE rscalar ;
    int scan = rb_scan_args(argc, argv, "11", &rb, &rscalar) ;

    if( scan == 1 )
    {
        if( rb_block_given_p() )
        {
            int i ;
            int size = num2size(rb) ;
            a = DMatrix_new(size, size) ;
            ra = wrap_dmatrix(a) ;
            
            VALUE result ;
        
            DMatrix_fill(a, dmatrix_doublereal_0) ;
            
            for( i = 0 ; i != size ; ++i )
            {
                result = rb_yield(INT2FIX(i)) ;
                a->data[i + i*a->vsize] = NUM2DBL(result) ;
            }
        }
        else if( rb_obj_is_kind_of(rb, rb_cDMatrix) )
        {
            int i ;
            DMatrix* b ;

            get_dmatrix(b, rb) ;
            
            if( b->hsize != 1 )
            {
                raise_dim_error() ;
            }
            
            a = DMatrix_new(b->vsize, b->vsize) ;
            ra = wrap_dmatrix(a) ;
            
            DMatrix_fill(a, dmatrix_doublereal_0) ;
            
            for( i = 0 ; i != b->vsize ; ++i )
            {
                a->data[i + i*a->vsize] = b->data[i] ;
            }
        }
        else if( TYPE(rb) == T_ARRAY )
        {
            VALUE* c_arr ;
            long size ;
            int i ;
            
            c_arr = RARRAY_PTR(rb) ;
            size = RARRAY_LEN(rb) ;
            
            if( size < 1 )
            {
                raise_dim_error() ;
            }
            
            a = DMatrix_new(size, size) ;
            ra = wrap_dmatrix(a) ;
            
            DMatrix_fill(a, dmatrix_doublereal_0) ;
            
            for( i = 0 ; i != size ; ++i )
            {
                a->data[i + i*a->vsize] = NUM2DBL(c_arr[i]) ;
            }
        }
        else
        {
            rb_raise(rb_eRuntimeError, "not a vector or array") ;
        }
    }
    else
    {
        int size = num2size(rb) ;
        a = DMatrix_new(size, size) ;
        ra = wrap_dmatrix(a) ;
        DMatrix_scalar_diagonal(a, NUM2DBL(rscalar)) ;
    }

    return ra ;
}