/*
*+
* 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