|
#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 );
|
|
}
|
