You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
562 lines
16 KiB
562 lines
16 KiB
/**************************************************************************
|
|
* *
|
|
* Copyright (C) 1992, Silicon Graphics, Inc. *
|
|
* *
|
|
* These coded instructions, statements, and computer programs contain *
|
|
* unpublished proprietary information of Silicon Graphics, Inc., and *
|
|
* are protected by Federal copyright law. They may not be disclosed *
|
|
* to third parties or copied or duplicated in any form, in whole or *
|
|
* in part, without the prior written consent of Silicon Graphics, Inc. *
|
|
* *
|
|
**************************************************************************/
|
|
|
|
/*
|
|
* arctessellator.c++ - $Revision: 1.6 $
|
|
* Derrick Burns - 1991
|
|
*/
|
|
|
|
#include "glimport.h"
|
|
#include "mystdio.h"
|
|
#include "myassert.h"
|
|
#include "arctess.h"
|
|
#include "bufpool.h"
|
|
#include "simplema.h"
|
|
#include "bezierar.h"
|
|
#include "trimvert.h"
|
|
#include "trimpool.h"
|
|
|
|
#define NOELIMINATION
|
|
|
|
#define steps_function(large, small, rate) (max(1, 1+ (int) ((large-small)/rate)));
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* ArcTessellator - construct an ArcTessellator
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
ArcTessellator::ArcTessellator( TrimVertexPool& t, Pool& p )
|
|
: trimvertexpool(t), pwlarcpool(p)
|
|
{
|
|
}
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* ~ArcTessellator - destroy an ArcTessellator
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
ArcTessellator::~ArcTessellator( void )
|
|
{
|
|
}
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* bezier - construct a bezier arc and attach it to an Arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::bezier( Arc *arc, REAL s1, REAL s2, REAL t1, REAL t2 )
|
|
{
|
|
assert( arc != 0 );
|
|
assert( ! arc->isTessellated() );
|
|
|
|
#ifndef NDEBUG
|
|
switch( arc->getside() ) {
|
|
case arc_left:
|
|
assert( s1 == s2 );
|
|
assert( t2 < t1 );
|
|
break;
|
|
case arc_right:
|
|
assert( s1 == s2 );
|
|
assert( t1 < t2 );
|
|
break;
|
|
case arc_top:
|
|
assert( t1 == t2 );
|
|
assert( s2 < s1 );
|
|
break;
|
|
case arc_bottom:
|
|
assert( t1 == t2 );
|
|
assert( s1 < s2 );
|
|
break;
|
|
case arc_none:
|
|
(void) abort();
|
|
break;
|
|
}
|
|
#endif
|
|
|
|
TrimVertex *p = trimvertexpool.get(2);
|
|
arc->pwlArc = new(pwlarcpool) PwlArc( 2, p );
|
|
p[0].param[0] = s1;
|
|
p[0].param[1] = t1;
|
|
p[1].param[0] = s2;
|
|
p[1].param[1] = t2;
|
|
assert( (s1 == s2) || (t1 == t2) );
|
|
arc->setbezier();
|
|
}
|
|
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* pwl_left - construct a left boundary pwl arc and attach it to an arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::pwl_left( Arc *arc, REAL s, REAL t1, REAL t2, REAL rate )
|
|
{
|
|
assert( t2 < t1 );
|
|
|
|
/* if(rate <= 0.06) rate = 0.06;*/
|
|
/* int nsteps = 1 + (int) ((t1 - t2) / rate ); */
|
|
int nsteps = steps_function(t1, t2, rate);
|
|
|
|
|
|
REAL stepsize = (t1 - t2) / (REAL) nsteps;
|
|
|
|
TrimVertex *newvert = trimvertexpool.get( nsteps+1 );
|
|
for( int i = nsteps; i > 0; i-- ) {
|
|
newvert[i].param[0] = s;
|
|
newvert[i].param[1] = t2;
|
|
t2 += stepsize;
|
|
}
|
|
newvert[i].param[0] = s;
|
|
newvert[i].param[1] = t1;
|
|
|
|
arc->makeSide( new(pwlarcpool) PwlArc( nsteps+1, newvert ), arc_left );
|
|
}
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* pwl_right - construct a right boundary pwl arc and attach it to an arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::pwl_right( Arc *arc, REAL s, REAL t1, REAL t2, REAL rate )
|
|
{
|
|
assert( t1 < t2 );
|
|
|
|
/* if(rate <= 0.06) rate = 0.06;*/
|
|
|
|
/* int nsteps = 1 + (int) ((t2 - t1) / rate ); */
|
|
int nsteps = steps_function(t2,t1,rate);
|
|
REAL stepsize = (t2 - t1) / (REAL) nsteps;
|
|
|
|
TrimVertex *newvert = trimvertexpool.get( nsteps+1 );
|
|
for( int i = 0; i < nsteps; i++ ) {
|
|
newvert[i].param[0] = s;
|
|
newvert[i].param[1] = t1;
|
|
t1 += stepsize;
|
|
}
|
|
newvert[i].param[0] = s;
|
|
newvert[i].param[1] = t2;
|
|
|
|
arc->makeSide( new(pwlarcpool) PwlArc( nsteps+1, newvert ), arc_right );
|
|
}
|
|
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* pwl_top - construct a top boundary pwl arc and attach it to an arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::pwl_top( Arc *arc, REAL t, REAL s1, REAL s2, REAL rate )
|
|
{
|
|
assert( s2 < s1 );
|
|
|
|
/* if(rate <= 0.06) rate = 0.06;*/
|
|
|
|
/* int nsteps = 1 + (int) ((s1 - s2) / rate ); */
|
|
int nsteps = steps_function(s1,s2,rate);
|
|
REAL stepsize = (s1 - s2) / (REAL) nsteps;
|
|
|
|
TrimVertex *newvert = trimvertexpool.get( nsteps+1 );
|
|
for( int i = nsteps; i > 0; i-- ) {
|
|
newvert[i].param[0] = s2;
|
|
newvert[i].param[1] = t;
|
|
s2 += stepsize;
|
|
}
|
|
newvert[i].param[0] = s1;
|
|
newvert[i].param[1] = t;
|
|
|
|
arc->makeSide( new(pwlarcpool) PwlArc( nsteps+1, newvert ), arc_top );
|
|
}
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* pwl_bottom - construct a bottom boundary pwl arc and attach it to an arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::pwl_bottom( Arc *arc, REAL t, REAL s1, REAL s2, REAL rate )
|
|
{
|
|
assert( s1 < s2 );
|
|
|
|
/* if(rate <= 0.06) rate = 0.06;*/
|
|
|
|
/* int nsteps = 1 + (int) ((s2 - s1) / rate ); */
|
|
int nsteps = steps_function(s2,s1,rate);
|
|
REAL stepsize = (s2 - s1) / (REAL) nsteps;
|
|
|
|
TrimVertex *newvert = trimvertexpool.get( nsteps+1 );
|
|
for( int i = 0; i < nsteps; i++ ) {
|
|
newvert[i].param[0] = s1;
|
|
newvert[i].param[1] = t;
|
|
s1 += stepsize;
|
|
}
|
|
newvert[i].param[0] = s2;
|
|
newvert[i].param[1] = t;
|
|
|
|
arc->makeSide( new(pwlarcpool) PwlArc( nsteps+1, newvert ), arc_bottom );
|
|
}
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* pwl - construct a pwl arc and attach it to an arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::pwl( Arc *arc, REAL s1, REAL s2, REAL t1, REAL t2, REAL rate )
|
|
{
|
|
|
|
/* if(rate <= 0.06) rate = 0.06;*/
|
|
|
|
int snsteps = 1 + (int) (abs(s2 - s1) / rate );
|
|
int tnsteps = 1 + (int) (abs(t2 - t1) / rate );
|
|
int nsteps = max(1,max( snsteps, tnsteps ));
|
|
|
|
REAL sstepsize = (s2 - s1) / (REAL) nsteps;
|
|
REAL tstepsize = (t2 - t1) / (REAL) nsteps;
|
|
TrimVertex *newvert = trimvertexpool.get( nsteps+1 );
|
|
for( long i = 0; i < nsteps; i++ ) {
|
|
newvert[i].param[0] = s1;
|
|
newvert[i].param[1] = t1;
|
|
s1 += sstepsize;
|
|
t1 += tstepsize;
|
|
}
|
|
newvert[i].param[0] = s2;
|
|
newvert[i].param[1] = t2;
|
|
|
|
/* arc->makeSide( new(pwlarcpool) PwlArc( nsteps+1, newvert ), arc_bottom ); */
|
|
arc->pwlArc = new(pwlarcpool) PwlArc( nsteps+1, newvert );
|
|
|
|
arc->clearbezier();
|
|
arc->clearside( );
|
|
}
|
|
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* tessellateLinear - constuct a linear pwl arc and attach it to an Arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::tessellateLinear( Arc *arc, REAL geo_stepsize, REAL arc_stepsize, int isrational )
|
|
{
|
|
assert( arc->pwlArc == NULL );
|
|
REAL s1, s2, t1, t2;
|
|
REAL stepsize = geo_stepsize * arc_stepsize;
|
|
BezierArc *b = arc->bezierArc;
|
|
|
|
if( isrational ) {
|
|
s1 = b->cpts[0] / b->cpts[2];
|
|
t1 = b->cpts[1] / b->cpts[2];
|
|
s2 = b->cpts[b->stride+0] / b->cpts[b->stride+2];
|
|
t2 = b->cpts[b->stride+1] / b->cpts[b->stride+2];
|
|
} else {
|
|
s1 = b->cpts[0];
|
|
t1 = b->cpts[1];
|
|
s2 = b->cpts[b->stride+0];
|
|
t2 = b->cpts[b->stride+1];
|
|
}
|
|
if( s1 == s2 )
|
|
if( t1 < t2 )
|
|
pwl_right( arc, s1, t1, t2, stepsize );
|
|
else
|
|
pwl_left( arc, s1, t1, t2, stepsize );
|
|
else if( t1 == t2 )
|
|
if( s1 < s2 )
|
|
pwl_bottom( arc, t1, s1, s2, stepsize );
|
|
else
|
|
pwl_top( arc, t1, s1, s2, stepsize );
|
|
else
|
|
pwl( arc, s1, s2, t1, t2, stepsize );
|
|
}
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* tessellateNonlinear - constuct a nonlinear pwl arc and attach it to an Arc
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
|
|
void
|
|
ArcTessellator::tessellateNonlinear( Arc *arc, REAL geo_stepsize, REAL arc_stepsize, int isrational )
|
|
{
|
|
assert( arc->pwlArc == NULL );
|
|
|
|
REAL stepsize = geo_stepsize * arc_stepsize;
|
|
|
|
int nsteps = 1 + (int) (1.0/stepsize);
|
|
|
|
TrimVertex *vert = trimvertexpool.get( nsteps+1 );
|
|
REAL dp = 1.0/nsteps;
|
|
BezierArc *bezierArc = arc->bezierArc;
|
|
|
|
arc->pwlArc = new(pwlarcpool) PwlArc();
|
|
arc->pwlArc->pts = vert;
|
|
|
|
if( isrational ) {
|
|
REAL pow_u[MAXORDER], pow_v[MAXORDER], pow_w[MAXORDER];
|
|
trim_power_coeffs( bezierArc, pow_u, 0 );
|
|
trim_power_coeffs( bezierArc, pow_v, 1 );
|
|
trim_power_coeffs( bezierArc, pow_w, 2 );
|
|
|
|
/* compute first point exactly */
|
|
REAL *b = bezierArc->cpts;
|
|
vert->param[0] = b[0]/b[2];
|
|
vert->param[1] = b[1]/b[2];
|
|
|
|
/* strength reduction on p = dp * step would introduce error */
|
|
int step;
|
|
int ocanremove = 0;
|
|
register long order = bezierArc->order;
|
|
for( step=1, ++vert; step<nsteps; step++, vert++ ) {
|
|
register REAL p = dp * step;
|
|
register REAL u = pow_u[0];
|
|
register REAL v = pow_v[0];
|
|
register REAL w = pow_w[0];
|
|
for( register int i = 1; i < order; i++ ) {
|
|
u = u * p + pow_u[i];
|
|
v = v * p + pow_v[i];
|
|
w = w * p + pow_w[i];
|
|
}
|
|
vert->param[0] = u/w;
|
|
vert->param[1] = v/w;
|
|
#ifndef NOELIMINATION
|
|
REAL ds = abs(vert[0].param[0] - vert[-1].param[0]);
|
|
REAL dt = abs(vert[0].param[1] - vert[-1].param[1]);
|
|
int canremove = (ds<geo_stepsize && dt<geo_stepsize) ? 1 : 0;
|
|
REAL ods=0.0, odt=0.0;
|
|
|
|
if( ocanremove && canremove ) {
|
|
REAL nds = ds + ods;
|
|
REAL ndt = dt + odt;
|
|
if( nds<geo_stepsize && ndt<geo_stepsize ) {
|
|
// remove previous point
|
|
--vert;
|
|
vert[0].param[0] = vert[1].param[0];
|
|
vert[0].param[1] = vert[1].param[1];
|
|
ods = nds;
|
|
odt = ndt;
|
|
ocanremove = 1;
|
|
} else {
|
|
ocanremove = canremove;
|
|
ods = ds;
|
|
odt = dt;
|
|
}
|
|
} else {
|
|
ocanremove = canremove;
|
|
ods = ds;
|
|
odt = dt;
|
|
}
|
|
#endif
|
|
}
|
|
|
|
/* compute last point exactly */
|
|
b += (order - 1) * bezierArc->stride;
|
|
vert->param[0] = b[0]/b[2];
|
|
vert->param[1] = b[1]/b[2];
|
|
|
|
} else {
|
|
REAL pow_u[MAXORDER], pow_v[MAXORDER];
|
|
trim_power_coeffs( bezierArc, pow_u, 0 );
|
|
trim_power_coeffs( bezierArc, pow_v, 1 );
|
|
|
|
/* compute first point exactly */
|
|
REAL *b = bezierArc->cpts;
|
|
vert->param[0] = b[0];
|
|
vert->param[1] = b[1];
|
|
|
|
/* strength reduction on p = dp * step would introduce error */
|
|
int step;
|
|
int ocanremove = 0;
|
|
register long order = bezierArc->order;
|
|
for( step=1, ++vert; step<nsteps; step++, vert++ ) {
|
|
register REAL p = dp * step;
|
|
register REAL u = pow_u[0];
|
|
register REAL v = pow_v[0];
|
|
for( register int i = 1; i < bezierArc->order; i++ ) {
|
|
u = u * p + pow_u[i];
|
|
v = v * p + pow_v[i];
|
|
}
|
|
vert->param[0] = u;
|
|
vert->param[1] = v;
|
|
#ifndef NOELIMINATION
|
|
REAL ds = abs(vert[0].param[0] - vert[-1].param[0]);
|
|
REAL dt = abs(vert[0].param[1] - vert[-1].param[1]);
|
|
int canremove = (ds<geo_stepsize && dt<geo_stepsize) ? 1 : 0;
|
|
REAL ods=0.0, odt=0.0;
|
|
|
|
if( ocanremove && canremove ) {
|
|
REAL nds = ds + ods;
|
|
REAL ndt = dt + odt;
|
|
if( nds<geo_stepsize && ndt<geo_stepsize ) {
|
|
// remove previous point
|
|
--vert;
|
|
vert[0].param[0] = vert[1].param[0];
|
|
vert[0].param[1] = vert[1].param[1];
|
|
ods = nds;
|
|
odt = ndt;
|
|
ocanremove = 1;
|
|
} else {
|
|
ocanremove = canremove;
|
|
ods = ds;
|
|
odt = dt;
|
|
}
|
|
} else {
|
|
ocanremove = canremove;
|
|
ods = ds;
|
|
odt = dt;
|
|
}
|
|
#endif
|
|
}
|
|
|
|
/* compute last point exactly */
|
|
b += (order - 1) * bezierArc->stride;
|
|
vert->param[0] = b[0];
|
|
vert->param[1] = b[1];
|
|
}
|
|
arc->pwlArc->npts = vert - arc->pwlArc->pts + 1;
|
|
/*
|
|
for( TrimVertex *vt=pwlArc->pts; vt != vert-1; vt++ ) {
|
|
if( tooclose( vt[0].param[0], vt[1].param[0] ) )
|
|
vt[1].param[0] = vt[0].param[0];
|
|
if( tooclose( vt[0].param[1], vt[1].param[1] ) )
|
|
vt[1].param[1] = vt[0].param[1];
|
|
}
|
|
*/
|
|
}
|
|
|
|
#ifdef NT
|
|
REAL ArcTessellator::gl_Bernstein[][MAXORDER][MAXORDER] = {
|
|
#else
|
|
const REAL ArcTessellator::gl_Bernstein[][MAXORDER][MAXORDER] = {
|
|
#endif
|
|
{
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{-1, 1, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{1, -2, 1, 0, 0, 0, 0, 0 },
|
|
{-2, 2, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{-1, 3, -3, 1, 0, 0, 0, 0 },
|
|
{3, -6, 3, 0, 0, 0, 0, 0 },
|
|
{-3, 3, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{1, -4, 6, -4, 1, 0, 0, 0 },
|
|
{-4, 12, -12, 4, 0, 0, 0, 0 },
|
|
{6, -12, 6, 0, 0, 0, 0, 0 },
|
|
{-4, 4, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{-1, 5, -10, 10, -5, 1, 0, 0 },
|
|
{5, -20, 30, -20, 5, 0, 0, 0 },
|
|
{-10, 30, -30, 10, 0, 0, 0, 0 },
|
|
{10, -20, 10, 0, 0, 0, 0, 0 },
|
|
{-5, 5, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{1, -6, 15, -20, 15, -6, 1, 0 },
|
|
{-6, 30, -60, 60, -30, 6, 0, 0 },
|
|
{15, -60, 90, -60, 15, 0, 0, 0 },
|
|
{-20, 60, -60, 20, 0, 0, 0, 0 },
|
|
{15, -30, 15, 0, 0, 0, 0, 0 },
|
|
{-6, 6, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 },
|
|
{0, 0, 0, 0, 0, 0, 0, 0 }
|
|
},
|
|
{
|
|
{-1, 7, -21, 35, -35, 21, -7, 1 },
|
|
{7, -42, 105, -140, 105, -42, 7, 0 },
|
|
{-21, 105, -210, 210, -105, 21, 0, 0 },
|
|
{35, -140, 210, -140, 35, 0, 0, 0 },
|
|
{-35, 105, -105, 35, 0, 0, 0, 0 },
|
|
{21, -42, 21, 0, 0, 0, 0, 0 },
|
|
{-7, 7, 0, 0, 0, 0, 0, 0 },
|
|
{1, 0, 0, 0, 0, 0, 0, 0 }
|
|
}};
|
|
|
|
|
|
/*-----------------------------------------------------------------------------
|
|
* trim_power_coeffs - compute power basis coefficients from bezier coeffients
|
|
*-----------------------------------------------------------------------------
|
|
*/
|
|
void
|
|
ArcTessellator::trim_power_coeffs( BezierArc *bez_arc, REAL *p, int coord )
|
|
{
|
|
register int stride = bez_arc->stride;
|
|
register int order = bez_arc->order;
|
|
register REAL *base = bez_arc->cpts + coord;
|
|
|
|
#ifdef NT
|
|
REAL (*mat)[MAXORDER][MAXORDER] = &gl_Bernstein[order-1];
|
|
REAL (*lrow)[MAXORDER] = &(*mat)[order];
|
|
|
|
REAL (*row)[MAXORDER] = &(*mat)[0];
|
|
|
|
for( ; row != lrow; row++ ) {
|
|
register REAL s = (REAL)0.0;
|
|
register REAL *point = base;
|
|
register REAL *mlast = *row + order;
|
|
for( REAL *m = *row; m != mlast; m++, point += stride )
|
|
s += *(m) * (*point);
|
|
*(p++) = s;
|
|
}
|
|
#else
|
|
REAL const (*mat)[MAXORDER][MAXORDER] = &gl_Bernstein[order-1];
|
|
REAL const (*lrow)[MAXORDER] = &(*mat)[order];
|
|
|
|
for( REAL const (*row)[MAXORDER] = &(*mat)[0]; row != lrow; row++ ) {
|
|
register REAL s = 0.0;
|
|
register REAL *point = base;
|
|
register REAL const *mlast = *row + order;
|
|
for( REAL const *m = *row; m != mlast; m++, point += stride )
|
|
s += *(m) * (*point);
|
|
*(p++) = s;
|
|
}
|
|
#endif
|
|
}
|