## File: palDmat.c

package info (click to toggle)
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182` ``````/* *+ * Name: * palDmat * Purpose: * Matrix inversion & solution of simultaneous equations * Language: * Starlink ANSI C * Type of Module: * Library routine * Invocation: * void palDmat( int n, double *a, double *y, double *d, int *jf, * int *iw ); * Arguments: * n = int (Given) * Number of simultaneous equations and number of unknowns. * a = double[] (Given & Returned) * A non-singular NxN matrix (implemented as a contiguous block * of memory). After calling this routine "a" contains the * inverse of the matrix. * y = double[] (Given & Returned) * On input the vector of N knowns. On exit this vector contains the * N solutions. * d = double * (Returned) * The determinant. * jf = int * (Returned) * The singularity flag. If the matrix is non-singular, jf=0 * is returned. If the matrix is singular, jf=-1 & d=0.0 are * returned. In the latter case, the contents of array "a" on * return are undefined. * iw = int[] (Given) * Integer workspace of size N. * Description: * Matrix inversion & solution of simultaneous equations * For the set of n simultaneous equations in n unknowns: * A.Y = X * this routine calculates the inverse of A, the determinant * of matrix A and the vector of N unknowns. * Authors: * PTW: Pat Wallace (STFC) * TIMJ: Tim Jenness (JAC, Hawaii) * {enter_new_authors_here} * History: * 2012-02-11 (TIMJ): * Combination of a port of the Fortran and a comparison * with the obfuscated GPL C routine. * Adapted with permission from the Fortran SLALIB library. * {enter_further_changes_here} * Notes: * - Implemented using Gaussian elimination with partial pivoting. * - Optimized for speed rather than accuracy with errors 1 to 4 * times those of routines optimized for accuracy. * Copyright: * Copyright (C) 2001 Rutherford Appleton Laboratory. * Copyright (C) 2012 Science and Technology Facilities Council. * All Rights Reserved. * Licence: * This program is free software: you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation, either * version 3 of the License, or (at your option) 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 Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General * License along with this program. If not, see * . * Bugs: * {note_any_bugs_here} *- */ #include "pal.h" void palDmat ( int n, double *a, double *y, double *d, int *jf, int *iw ) { const double SFA = 1e-20; int k; double*aoff; *jf=0; *d=1.0; for(k=0,aoff=a; kamx){ amx=t; imx=i; aoff2=apos2; } } } if(amx0;){ int ki=iw[k]; if(k!=ki){ int i; double *apos = a; for(i=0;i