home *** CD-ROM | disk | FTP | other *** search
- // VirtualDub - Video processing and capture application
- // System library component
- // Copyright (C) 1998-2004 Avery Lee, All Rights Reserved.
- //
- // Beginning with 1.6.0, the VirtualDub system library is licensed
- // differently than the remainder of VirtualDub. This particular file is
- // thus licensed as follows (the "zlib" license):
- //
- // This software is provided 'as-is', without any express or implied
- // warranty. In no event will the authors be held liable for any
- // damages arising from the use of this software.
- //
- // Permission is granted to anyone to use this software for any purpose,
- // including commercial applications, and to alter it and redistribute it
- // freely, subject to the following restrictions:
- //
- // 1. The origin of this software must not be misrepresented; you must
- // not claim that you wrote the original software. If you use this
- // software in a product, an acknowledgment in the product
- // documentation would be appreciated but is not required.
- // 2. Altered source versions must be plainly marked as such, and must
- // not be misrepresented as being the original software.
- // 3. This notice may not be removed or altered from any source
- // distribution.
-
- #include "stdafx.h"
- #include <vd2/system/vdstl.h>
- #include <vd2/system/vectors.h>
-
- bool VDSolveLinearEquation(double *src, int n, ptrdiff_t stride_elements, double *b, double tolerance) {
- vdfastvector<double *> array(n);
- double **m = &array[0];
- int i, j, k;
-
- for(i=0; i<n; ++i) {
- m[i] = src;
- src += stride_elements;
- }
-
- // factor U
- for(i=0; i<n; ++i) {
- int best = i;
-
- for(j=i+1; j<n; ++j) {
- if (fabs(m[best][i]) < fabs(m[j][i]))
- best = j;
- }
-
- std::swap(m[i], m[best]);
- std::swap(b[i], b[best]);
-
- if (fabs(m[i][i]) < tolerance)
- return false;
-
- double f = 1.0 / m[i][i];
-
- m[i][i] = 1.0;
-
- for(j=i+1; j<n; ++j)
- m[i][j] *= f;
-
- b[i] *= f;
-
- for(j=i+1; j<n; ++j) {
- b[j] -= b[i] * m[j][i];
- for(k=n-1; k>=i; --k)
- m[j][k] -= m[i][k] * m[j][i];
- }
- }
-
- // factor L
- for(i=n-1; i>=0; --i)
- for(j=i-1; j>=0; --j)
- b[j] -= b[i] * m[j][i];
-
- return true;
- }
-