[prev in list] [next in list] [prev in thread] [next in thread]
List: koffice-devel
Subject: Linear algebra library
From: Inge Wallin <inge () lysator ! liu ! se>
Date: 2006-07-27 14:14:31
Message-ID: 200607271614.31615.inge () lysator ! liu ! se
[Download RAW message or body]
On Wednesday 26 July 2006 20.41, Benoît Jacob wrote:
> Hi,
>
> after receiving many suggestions, especially in comments on Carsten's blog,
> and also googling a lot, I evaluated several linear algebra packages, and
> my feeling is that we should use GSL, the GNU Scientific Library:
>
> http://www.gnu.org/software/gsl/
After a length discussion on IRC, where size considerations were raised over
GSL, ublas, and other comprehensive, but big libraries, it was decided to
look into actually including a small library into the apps. An alternative
would be to try to get it into kdelibs or even Qt, but that is a longer term
venture.
Naturally it's madness to write a small matrix / linalg library from scratch
since that has been done 100s of times. So I dug out some old code of mine
and a friends that was used in a graphics library called SiPP, the Simple
Polygon Processor.
Here it is, it's C not C++, it is funnily indented, but it's very well tested,
it's fast and it's readable. It should be easy both to C++-ify and to extend
with some needed functions. Do with it what you want.
-Inge
> it looks very well-done, provides all what we might possibly need (not only
> linear algebra, but a lot more math stuff), is widely used, and is
> portable: http://www.gnu.org/software/gsl/#platforms
>
> Under Windows, is it OK to use minGW or Cygwin ? Or is it necessary to only
> use libs that can be compiled with MSVC++? In that case, according to this
> link:
>
> http://www.mail-archive.com/bug-gsl%40gnu.org/msg00151.html
>
> it should be possible to build with MSVC++ .Net or later versions, but one
> might be out of luck with MSVC++ 6.0 (not sure).
>
> Other options included Boost/uBLAS, but unfortunately it doesn't have
> enough features so we'd end up reinventing the wheel very often. For
> instance it can't compute eigenvectors. Another option was TNT/JAMA++, but
> its website wasn't very comforting, with the source code only available as
> a .zip archive... Other small projects like Newmat looked good but we're
> looking for a big well-tested project aren't we, if we ever want to propose
> to OpenBabel to use it.
>
> The only drawback of GSL is that it's C-only. Fortunately, there are some
> C++ wrappers for GSL out there. The one that looks most promising and
> actively developed is GSL--:
>
> http://cholm.home.cern.ch/cholm/misc/#gslmm
>
> I emailed its author, who told me it's currently 25% complete, adding that
> he'd be glad to accept patches. I think that such patches should be easy to
> write as it's only wrapper stuff, the functions themselves are already
> provided by GSL.
>
> So here's what I propose to do:
>
> 0) I'll be away from July 30 to August 12, so I'll only start then.
> 1) Begin using GSL functions in Kalzium where needed, for instance to
> compute the plane best approximating n atoms.
> 2) Port that to GSL--, sending patches to GSL-- as required.
> 3) Completely port Kalzium to GSL--, again sending patches to GSL-- as
> required. That means for instance that instead of using OpenBabel's vector3
> class, we'd now use its equivalent in GSL--.
> 4) Propose to OpenBabel patches turning the vector3, matrix3x3, ... classes
> into wrappers around GSL-- stuff. Using C++ inheritance, we could do that
> while preserving compatibility. For instance, momentarily call GSL::vector3
> the class in GSL-- that represents a 3d vector of doubles. Then we could
> replace OpenBabel::vector3 by
>
> class OpenBabel::vector3 : public GSL::vector3
> {
> // here, add the methods that are in the current OpenBabel::vector3
> // class, but that are not in GSL::vector3, like createOrthoVector
>
> // for many methods, like the operators, there should be nothing to
> // do, though
> }
>
> What do you think? This way Kalzium and/or OpenBabel wouldn't have to worry
> anymore about math stuff, risk of floating-point violations, etc...
>
> Cheers,
> Benoit
> _______________________________________________
> Kalzium mailing list
> Kalzium@kde.org
> https://mail.kde.org/mailman/listinfo/kalzium
--
Inge Wallin | Thus spake the master programmer: |
| "After three days without programming, |
inge@lysator.liu.se | life becomes meaningless." |
| Geoffrey James: The Tao of Programming. |
["geometric.c" (text/x-csrc)]
/**
** sipp - SImple Polygon Processor
**
** A general 3d graphic package
**
** Copyright Equivalent Software HB 1992
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation; either version 1, or any later version.
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
** You can receive a copy of the GNU General Public License from the
** Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**/
/**
** geometric.c - Matrixes, transformations and coordinates.
**/
#include <stdio.h>
#include <math.h>
#include <sipp.h>
#include <smalloc.h>
#include <geometric.h>
/* =================================================================== */
/* */
Transf_mat ident_matrix = {{ /* Unit tranfs. matrix */
{ 1.0, 0.0, 0.0 },
{ 0.0, 1.0, 0.0 },
{ 0.0, 0.0, 1.0 },
{ 0.0, 0.0, 0.0 }
}};
/* =================================================================== */
/*
* Allocate a new matrix, and if INITMAT != NULL copy the contents
* of INITMAT to the new matrix, otherwise copy the identity matrix
* to the new matrix.
*/
Transf_mat *
transf_mat_create(initmat)
Transf_mat * initmat;
{
Transf_mat * mat;
mat = (Transf_mat *) smalloc(sizeof(Transf_mat));
if (initmat != NULL)
MatCopy(mat, initmat);
else
MatCopy(mat, &ident_matrix);
return mat;
}
void
transf_mat_destruct(mat)
Transf_mat * mat;
{
sfree(mat);
}
/* =================================================================== */
/* Transformation routines (see also geometric.h) */
/*
* Normalize a vector.
*/
void
vecnorm(vec)
Vector *vec;
{
double len;
len = VecLen(*vec);
if (len == 0.0) { /* As Mark says, we could really use error handling...*/
MakeVector(*vec, 0.0, 0.0, 0.0);
} else {
VecScalMul(*vec, 1.0 / len, *vec);
}
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a translation along the vector described by DX, DY and DZ.
*
* [a b c 0] [ 1 0 0 0] [ a b c 0]
* [d e f 0] [ 0 1 0 0] [ d e f 0]
* [g h i 0] * [ 0 0 1 0] = [ g h i 0]
* [j k l 1] [Tx Ty Tz 1] [j+Tx k+Ty l+Tz 1]
*/
void
mat_translate(mat, dx, dy, dz)
Transf_mat * mat;
double dx;
double dy;
double dz;
{
mat->mat[3][0] += dx;
mat->mat[3][1] += dy;
mat->mat[3][2] += dz;
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a rotation with the angle ANG around the X axis.
*
* [a b c 0] [1 0 0 0] [a b*Ca-c*Sa b*Sa+c*Ca 0]
* [d e f 0] [0 Ca Sa 0] [d e*Ca-f*Sa e*Sa+f*Ca 0]
* [g h i 0] * [0 -Sa Ca 0] = [g h*Ca-i*Sa h*Sa+i*Ca 0]
* [j k l 1] [0 0 0 1] [j k*Ca-l*Sa k*Se+l*Ca 1]
*/
void
mat_rotate_x(mat, ang)
Transf_mat * mat;
double ang;
{
double cosang;
double sinang;
double tmp;
int i;
cosang = cos(ang);
sinang = sin(ang);
if (fabs(cosang) < 1.0e-15) {
cosang = 0.0;
}
if (fabs(sinang) < 1.0e-15) {
sinang = 0.0;
}
for (i = 0; i < 4; ++i) {
tmp = mat->mat[i][1];
mat->mat[i][1] = mat->mat[i][1] * cosang
- mat->mat[i][2] * sinang;
mat->mat[i][2] = tmp * sinang + mat->mat[i][2] * cosang;
}
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a rotation with the angle ANG around the Y axis.
*
* [a b c 0] [Ca 0 -Sa 0] [a*Ca+c*Sa b -a*Sa+c*Ca 0]
* [d e f 0] * [ 0 1 0 0] = [d*Ca+f*Sa e -d*Sa+f*Ca 0]
* [g h i 0] [Sa 0 Ca 0] [g*Ca+i*Sa h -g*Sa+i*Ca 0]
* [j k l 1] [ 0 0 0 1] [j*Ca+l*Sa k -j*Sa+l*Ca 1]
*/
void
mat_rotate_y(mat, ang)
Transf_mat * mat;
double ang;
{
double cosang;
double sinang;
double tmp;
int i;
cosang = cos(ang);
sinang = sin(ang);
if (fabs(cosang) < 1.0e-15) {
cosang = 0.0;
}
if (fabs(sinang) < 1.0e-15) {
sinang = 0.0;
}
for (i = 0; i < 4; ++i) {
tmp = mat->mat[i][0];
mat->mat[i][0] = mat->mat[i][0] * cosang
+ mat->mat[i][2] * sinang;
mat->mat[i][2] = -tmp * sinang + mat->mat[i][2] * cosang;
}
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a rotation with the angle ANG around the Z axis.
*
* [a b c 0] [ Ca Sa 0 0] [a*Ca-b*Sa a*Sa+b*Ca c 0]
* [d e f 0] [-Sa Ca 0 0] [d*Ca-e*Sa d*Sa+e*Ca f 0]
* [g h i 0] * [ 0 0 1 0] = [g*Ca-h*Sa g*Sa+h*Ca i 0]
* [j k l 1] [ 0 0 0 1] [j*Ca-k*Sa j*Sa+k*Ca l 0]
*/
void
mat_rotate_z(mat, ang)
Transf_mat * mat;
double ang;
{
double cosang;
double sinang;
double tmp;
int i;
cosang = cos(ang);
sinang = sin(ang);
if (fabs(cosang) < 1.0e-15) {
cosang = 0.0;
}
if (fabs(sinang) < 1.0e-15) {
sinang = 0.0;
}
for (i = 0; i < 4; ++i) {
tmp = mat->mat[i][0];
mat->mat[i][0] = mat->mat[i][0] * cosang
- mat->mat[i][1] * sinang;
mat->mat[i][1] = tmp * sinang + mat->mat[i][1] * cosang;
}
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a rotation with the angle ANG around the line represented
* by the point POINT and the vector VECTOR.
*/
void
mat_rotate(mat, point, vector, ang)
Transf_mat * mat;
Vector * point;
Vector * vector;
double ang;
{
double ang2;
double ang3;
ang2 = atan2(vector->y, vector->x);
ang3 = atan2(hypot(vector->x, vector->y), vector->z);
mat_translate(mat, -point->x, -point->y, -point->z);
mat_rotate_z(mat, -ang2);
mat_rotate_y(mat, -ang3);
mat_rotate_z(mat, ang);
mat_rotate_y(mat, ang3);
mat_rotate_z(mat, ang2);
mat_translate(mat, point->x, point->y, point->z);
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a scaling with the scaling factors XSCALE, YSCALE and ZSCALE
*
* [a b c 0] [Sx 0 0 0] [a*Sx b*Sy c*Sz 0]
* [d e f 0] [ 0 Sy 0 0] [d*Sx e*Sy f*Sz 0]
* [g h i 0] * [ 0 0 Sz 0] = [g*Sx h*Sy i*Sz 0]
* [j k l 1] [ 0 0 0 1] [j*Sx k*Sy l*Sz 1]
*/
void
mat_scale(mat, xscale, yscale, zscale)
Transf_mat * mat;
double xscale;
double yscale;
double zscale;
{
int i;
for (i = 0; i < 4; ++i) {
mat->mat[i][0] *= xscale;
mat->mat[i][1] *= yscale;
mat->mat[i][2] *= zscale;
}
}
/*
* Set MAT to the transformation matrix that represents the
* concatenation between the previous transformation in MAT
* and a mirroring in the plane defined by the point POINT
* and the normal vector NORM.
*/
void
mat_mirror_plane(mat, point, norm)
Transf_mat * mat;
Vector * point;
Vector * norm;
{
Transf_mat tmp;
double factor;
/* The first thing we do is to make a transformation matrix */
/* for mirroring through a plane with the same normal vector */
/* as our, but through the origin instead. */
factor = 2.0 / (norm->x * norm->x + norm->y * norm->y
+ norm->z * norm->z);
/* The diagonal elements. */
tmp.mat[0][0] = 1 - factor * norm->x * norm->x;
tmp.mat[1][1] = 1 - factor * norm->y * norm->y;
tmp.mat[2][2] = 1 - factor * norm->z * norm->z;
/* The rest of the matrix */
tmp.mat[1][0] = tmp.mat[0][1] = -factor * norm->x * norm->y;
tmp.mat[2][0] = tmp.mat[0][2] = -factor * norm->x * norm->z;
tmp.mat[2][1] = tmp.mat[1][2] = -factor * norm->y * norm->z;
tmp.mat[3][0] = tmp.mat[3][1] = tmp.mat[3][2] = 0.0;
/* Do the actual transformation. This is done in 3 steps: */
/* 1) Translate the plane so that it goes through the origin. */
/* 2) Do the actual mirroring. */
/* 3) Translate it all back to the starting position. */
mat_translate(mat, -point->x, -point->y, -point->z);
mat_mul(mat, mat, &tmp);
mat_translate(mat, point->x, point->y, point->z);
}
/*
* Multiply the Matrix A with the Matrix B, and store the result
* into the Matrix RES. It is possible for RES to point to the
* same Matrix as either A or B since the result is stored into
* a temporary during computation.
*
* [a b c 0] [A B C 0] [aA+bD+cG aB+bE+cH aC+bF+cI 0]
* [d e f 0] [D E F 0] [dA+eD+fG dB+eE+fH dC+eF+fI 0]
* [g h i 0] [G H I 0] = [gA+hD+iG gB+hE+iH gC+hF+iI 0]
* [j k l 1] [J K L 1] [jA+kD+lG+J jB+kE+lH+K jC+kF+lI+L 1]
*/
void
mat_mul(res, a, b)
Transf_mat * res;
Transf_mat * a;
Transf_mat * b;
{
Transf_mat tmp;
int i;
for (i = 0; i < 4; ++i) {
tmp.mat[i][0] = a->mat[i][0] * b->mat[0][0]
+ a->mat[i][1] * b->mat[1][0]
+ a->mat[i][2] * b->mat[2][0];
tmp.mat[i][1] = a->mat[i][0] * b->mat[0][1]
+ a->mat[i][1] * b->mat[1][1]
+ a->mat[i][2] * b->mat[2][1];
tmp.mat[i][2] = a->mat[i][0] * b->mat[0][2]
+ a->mat[i][1] * b->mat[1][2]
+ a->mat[i][2] * b->mat[2][2];
}
tmp.mat[3][0] += b->mat[3][0];
tmp.mat[3][1] += b->mat[3][1];
tmp.mat[3][2] += b->mat[3][2];
MatCopy(res, &tmp);
}
/*
* Transform the Point3d VEC with the transformation matrix MAT, and
* put the result into the vector *RES.
*
* [a b c 0]
* [d e f 0]
* [x y z 1] [g h i 0] = [ax+dy+gz+j bx+ey+hz+k cx+fy+iz+l 1]
* [j k l 1]
*/
void
point_transform(res, vec, mat)
Vector * res;
Vector * vec;
Transf_mat * mat;
{
res->x = mat->mat[0][0] * vec->x + mat->mat[1][0] * vec->y
+ mat->mat[2][0] * vec->z + mat->mat[3][0];
res->y = mat->mat[0][1] * vec->x + mat->mat[1][1] * vec->y
+ mat->mat[2][1] * vec->z + mat->mat[3][1];
res->z = mat->mat[0][2] * vec->x + mat->mat[1][2] * vec->y
+ mat->mat[2][2] * vec->z + mat->mat[3][2];
}
["geometric.h" (text/x-chdr)]
/**
** sipp - SImple Polygon Processor
**
** A general 3d graphic package
**
** Copyright Equivalent Software HB 1992
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation; either version 1, or any later version.
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
** You can receive a copy of the GNU General Public License from the
** Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**/
/**
** geometric.h - All kinds of stuff with matrixes, transformations,
** coordinates
**/
/* Make sure no multiple including */
#ifndef _GEOMETRIC_H_
#define _GEOMETRIC_H_
#include <math.h>
#ifndef NOMEMCPY
#include <memory.h>
#endif
/* #define PI 3.1415926535897932384626 */
typedef struct {
double x, y, z;
} Vector;
/*
* NOTE:
* Capitalized types denote Vectors and other aggregates.
* Lower case denote scalars.
*/
/* V = Vec(x, y, z) */
#define MakeVector(V, xx, yy, zz) { (V).x=(xx); \
(V).y=(yy); \
(V).z=(zz); }
/* A = -A */
#define VecNegate(A) { (A).x=0-(A).x; \
(A).y=0-(A).y; \
(A).z=0-(A).z; }
/* return A . B */
#define VecDot(A, B) ((A).x*(B).x+(A).y*(B).y+(A).z*(B).z)
/* return length(A) */
#define VecLen(A) (sqrt((double)VecDot(A, A)))
/* B = A */
#define VecCopy(B, A) ((B) = (A))
/* C = A + B */
#define VecAdd(C, A, B) { (C).x=(A).x+(B).x; \
(C).y=(A).y+(B).y; \
(C).z=(A).z+(B).z; }
/* C = A - B */
#define VecSub(C, A, B) { (C).x=(A).x-(B).x; \
(C).y=(A).y-(B).y; \
(C).z=(A).z-(B).z; }
/* C = a*A */
#define VecScalMul(C, a, A) { (C).x=(a)*(A).x; \
(C).y=(a)*(A).y; \
(C).z=(a)*(A).z; }
/* C = a*A + B */
#define VecAddS(C, a, A, B) { (C).x=(a)*(A).x+(B).x; \
(C).y=(a)*(A).y+(B).y; \
(C).z=(a)*(A).z+(B).z; }
/* C = a*A + b*B */
#define VecComb(C, a, A, b, B) { (C).x=(a)*(A).x+(b)*(B).x; \
(C).y=(a)*(A).y+(b)*(B).y; \
(C).z=(a)*(A).z+(b)*(B).z; }
/* C = A X B */
#define VecCross(C, A, B) { (C).x=(A).y*(B).z-(A).z*(B).y; \
(C).y=(A).z*(B).x-(A).x*(B).z; \
(C).z=(A).x*(B).y-(A).y*(B).x; }
#define VecMax(C, A, B) { (C).x=(((A).x>(B).x)?(A).x:(B).x); \
(C).y=(((A).y>(B).y)?(A).y:(B).y); \
(C).z=(((A).z>(B).z)?(A).z:(B).z); }
#define VecMin(C, A, B) { (C).x=(((A).x<(B).x)?(A).x:(B).x); \
(C).y=(((A).y<(B).y)?(A).y:(B).y); \
(C).z=(((A).z<(B).z)?(A).z:(B).z); }
/* ================================================================ */
/* Matrix operations */
/*
* Define a homogenous transformation matrix. The first row (vector)
* is the new X axis, i.e. the X axis in the transformed coordinate
* system. The second row is the new Y axis, and so on. The last row
* is the translation, for a transformed point.
*
* The reason we make surround the rows with a struct is that we
* don't want to say (Transf_mat *) &foo[0] instead of &foo when
* sending an address to a matrix as a parameter to a function.
* Alas, arrays are not first class objects in C.
*/
typedef struct {
double mat[4][3];
} Transf_mat;
extern Transf_mat ident_matrix;
/* *A = *B N.b. A and B are pointers! */
#define MatCopy(A, B) (*A) = (*B)
/*----------------------------------------------------------------------*/
/* Function declarations for the functions in geometric.c */
EXTERN void
vecnorm _ANSI_ARGS_((Vector *vec));
EXTERN Transf_mat *
transf_mat_create _ANSI_ARGS_((Transf_mat *initmat));
EXTERN void
transf_mat_destruct _ANSI_ARGS_((Transf_mat *mat));
EXTERN void
mat_translate _ANSI_ARGS_((Transf_mat *mat,
double dx,
double dy,
double dz));
EXTERN void
mat_rotate_x _ANSI_ARGS_((Transf_mat *mat,
double ang));
EXTERN void
mat_rotate_y _ANSI_ARGS_((Transf_mat *mat,
double ang));
EXTERN void
mat_rotate_z _ANSI_ARGS_((Transf_mat *mat,
double ang));
EXTERN void
mat_rotate _ANSI_ARGS_((Transf_mat *mat,
Vector *point,
Vector *vector,
double ang));
EXTERN void
mat_scale _ANSI_ARGS_((Transf_mat *mat,
double xscale,
double yscale,
double zscale));
EXTERN void
mat_mirror_plane _ANSI_ARGS_((Transf_mat *mat,
Vector *point,
Vector *norm));
EXTERN void
mat_mul _ANSI_ARGS_((Transf_mat *res,
Transf_mat *a,
Transf_mat *b));
EXTERN void
point_transform _ANSI_ARGS_((Vector *res,
Vector *vec,
Transf_mat *mat));
#endif /* _GEOMETRIC_H_ */
_______________________________________________
koffice-devel mailing list
koffice-devel@kde.org
https://mail.kde.org/mailman/listinfo/koffice-devel
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic