SGI SI GLU library

1 files changed, 609 insertions, 0 deletions

diff --git a/src/glu/sgi/libnurbs/internals/arctess.cc b/src/glu/sgi/libnurbs/internals/arctess.cc
new file mode 100644
index 00000000000..0559c35029e
--- /dev/null
+++ b/src/glu/sgi/libnurbs/internals/arctess.cc
@@ -0,0 +1,609 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+**
+** http://oss.sgi.com/projects/FreeB
+**
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+**
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+**
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+*/
+
+/*
+ * arctessellator.c++
+ *
+ * $Date: 2001/03/17 00:25:40 $ $Revision: 1.1 $
+ * $Header: /home/krh/git/sync/mesa-cvs-repo/Mesa/src/glu/sgi/libnurbs/internals/arctess.cc,v 1.1 2001/03/17 00:25:40 brianp Exp $
+ */
+
+#include "glimports.h"
+#include "mystdio.h"
+#include "myassert.h"
+#include "arctess.h"
+#include "bufpool.h"
+#include "simplemath.h"
+#include "bezierarc.h"
+#include "trimvertex.h"
+#include "trimvertpool.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 );
+    int i;
+    for( 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 );
+    int i;
+    for( 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 );
+    int i;
+    for( 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 );
+    int i;
+    for( 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 );
+    long i;
+    for( 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;
+
+    //we don't need to scale by arc_stepsize if the trim curve
+    //is piecewise linear. Reason: In pwl_right, pwl_left, pwl_top, pwl_left,
+    //and pwl, the nsteps is computed by deltaU (or V) /stepsize. 
+    //The quantity deltaU/arc_stepsize doesn't have any meaning. And
+    //it causes problems: see bug 517641
+    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;
+
+    BezierArc *bezierArc = arc->bezierArc;
+
+    REAL size; //bounding box size of the curve in UV 
+    {
+      int i,j;
+      REAL min_u, min_v, max_u,max_v;
+      min_u = max_u = bezierArc->cpts[0];
+      min_v = max_v = bezierArc->cpts[1];
+      for(i=1, j=2; i<bezierArc->order; i++, j+= bezierArc->stride)
+	{
+	  if(bezierArc->cpts[j] < min_u)
+	    min_u = bezierArc->cpts[j];
+	  if(bezierArc->cpts[j] > max_u)
+	    max_u = bezierArc->cpts[j];
+	  if(bezierArc->cpts[j+1] < min_v)
+	    min_v = bezierArc->cpts[j+1];	  
+	  if(bezierArc->cpts[j+1] > max_v)
+	    max_v = bezierArc->cpts[j+1]; 
+	}
+
+      size = max_u - min_u;
+      if(size < max_v - min_v)
+	size = max_v - min_v;
+    }
+      
+    /*int	nsteps 		= 1 + (int) (1.0/stepsize);*/
+
+    int nsteps = (int) (size/stepsize);
+    if(nsteps <=0)
+      nsteps=1;
+
+    TrimVertex *vert	= trimvertexpool.get( nsteps+1 );
+    REAL dp 		= 1.0/nsteps;
+
+
+    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];
+    }
+*/
+}
+
+const REAL ArcTessellator::gl_Bernstein[][MAXORDER][MAXORDER] = {
+ {
+  {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;
+
+    REAL const (*mat)[MAXORDER][MAXORDER] = &gl_Bernstein[order-1];
+    REAL const (*lrow)[MAXORDER] = &(*mat)[order];
+
+    /* WIN32 didn't like the following line within the for-loop */
+    REAL const (*row)[MAXORDER] =  &(*mat)[0];
+    for( ; 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;
+    }
+}