From koffice-devel Thu Jul 27 14:14:31 2006 From: Inge Wallin Date: Thu, 27 Jul 2006 14:14:31 +0000 To: koffice-devel Subject: Linear algebra library Message-Id: <200607271614.31615.inge () lysator ! liu ! se> X-MARC-Message: https://marc.info/?l=koffice-devel&m=115401164714171 MIME-Version: 1 Content-Type: multipart/mixed; boundary="--Boundary-00=_HpMyEpz1r4XadF5" --Boundary-00=_HpMyEpz1r4XadF5 Content-Type: text/plain; charset="iso-8859-1" Content-Transfer-Encoding: quoted-printable Content-Disposition: inline On Wednesday 26 July 2006 20.41, Beno=EEt Jacob wrote: > Hi, > > after receiving many suggestions, especially in comments on Carsten's blo= g, > 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 ove= r=20 GSL, ublas, and other comprehensive, but big libraries, it was decided to=20 look into actually including a small library into the apps. An alternative= =20 would be to try to get it into kdelibs or even Qt, but that is a longer ter= m=20 venture. Naturally it's madness to write a small matrix / linalg library from scratc= h=20 since that has been done 100s of times. So I dug out some old code of mine= =20 and a friends that was used in a graphics library called SiPP, the Simple=20 Polygon Processor. Here it is, it's C not C++, it is funnily indented, but it's very well test= ed,=20 it's fast and it's readable. It should be easy both to C++-ify and to exten= d=20 with some needed functions. Do with it what you want. -Inge > it looks very well-done, provides all what we might possibly need (not on= ly > 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 on= ly > 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 propo= se > 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 vecto= r3 > class, we'd now use its equivalent in GSL--. > 4) Propose to OpenBabel patches turning the vector3, matrix3x3, ... class= es > into wrappers around GSL-- stuff. Using C++ inheritance, we could do that > while preserving compatibility. For instance, momentarily call GSL::vecto= r3 > 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 wor= ry > 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 =2D-=20 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. = | --Boundary-00=_HpMyEpz1r4XadF5 Content-Type: text/x-csrc; charset="iso-8859-1"; name="geometric.c" Content-Transfer-Encoding: 7bit Content-Disposition: attachment; filename="geometric.c" /** ** 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 #include #include #include #include /* =================================================================== */ /* */ 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]; } --Boundary-00=_HpMyEpz1r4XadF5 Content-Type: text/x-chdr; charset="iso-8859-1"; name="geometric.h" Content-Transfer-Encoding: 7bit Content-Disposition: attachment; filename="geometric.h" /** ** 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 #ifndef NOMEMCPY #include #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_ */ --Boundary-00=_HpMyEpz1r4XadF5 Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ koffice-devel mailing list koffice-devel@kde.org https://mail.kde.org/mailman/listinfo/koffice-devel --Boundary-00=_HpMyEpz1r4XadF5--