kv

-

kernel handle of type MKernelVector

rt

-

type ALGEB rtable object

i

-

ith indexing function in an rtable

    #include <string.h>

    #include "maplec.h"

    ALGEB M_DECL MySetDiagonal( MKernelVector kv, ALGEB *args )

    {

    M_INT argc, i, diag, size, index[2], num_iter;

    RTableSettings rts;

    RTableData rtd;

    ALGEB rt;

    argc = MapleNumArgs(kv,(ALGEB)args);

    if( argc != 3 ) {

        MapleRaiseError(kv,"three arguments expected");

            return( NULL );

        }

    rt = args[1];

    if( !IsMapleRTable(kv,rt) ) {

        MapleRaiseError(kv,"rtable expected");

            return( NULL );

        }

    RTableGetSettings(kv,&rts,rt);

    if( rts.data_type != RTABLE_FLOAT64 ) {

        MapleRaiseError(kv,"float[8] rtable expected");

            return( NULL );

        }

    if( rts.num_dimensions != 2 ) {

        MapleRaiseError(kv,"2D rtable expected");

            return( NULL );

        }

    if( RTableUpperBound(kv,rt,1) != RTableUpperBound(kv,rt,2)

        || RTableLowerBound(kv,rt,1) != 1

        || RTableLowerBound(kv,rt,2) != 1 ) {

        MapleRaiseError(kv,

                "square rtable with indexing starting at 1 expected");

            return( NULL );

        }

    rtd.float64 = MapleToFloat64(kv,args[3]);

    diag = MapleToM_INT(kv,args[2]);

    size = RTableUpperBound(kv,rt,1) - RTableLowerBound(kv,rt,1) + 1;

    if( diag > size-1 )

        return( rt );

    if( diag < 0 ) {

        index[0] = 1;

        index[1] = -diag+1;

        num_iter = size - index[1] + 1;

    }

    else {

        index[0] = diag+1;

        index[1] = 1;

        num_iter = size - index[0] + 1;

    }

    if( MapleNumArgs(kv,rts.index_functions) == 1

        && RTableIndFn(kv,rt,1) == RTABLE_INDEX_BAND

        && rts.order == RTABLE_FORTRAN

        && rts.storage == RTABLE_BAND ) {

        M_INT p1, p2, j;

        ALGEB band_args;

        FLOAT64 *data;

        band_args = RTableIndFnArgs(kv,rt,1);

        p1 = MapleToM_INT(kv,(ALGEB)band_args[1]);

        p2 = MapleToM_INT(kv,(ALGEB)band_args[2]);

        if( p1 == rts.p1 && p2 == rts.p2 ) {

        if( (diag < 0 && -diag > p2)

            || (diag > 0 && diag > p1) ) {

            MapleRaiseError(kv,"diagonal outside of bands");

                    return( NULL );

                }

        data = (FLOAT64*)RTableDataBlock(kv,rt);

        j = MATRIX_OFFSET_FORTRAN_BAND(index[0],index[1],

                           size,size,p1,p2);

        for( i=0; i<num_iter; ++i ) {

            data[j] = rtd.float64;

            j += p1+p2+1;

        }

        return( rt );

        }

    }

    for( i=0; i<num_iter; ++i ) {

        RTableAssign(kv,rt,index,rtd);

        index[0]++;

        index[1]++;

    }

    return( rt );

    }