mirror of
https://github.com/php/pecl-math-stats.git
synced 2026-03-23 23:22:18 +01:00
92 lines
2.6 KiB
C
92 lines
2.6 KiB
C
#include <math.h>
|
|
float sdot(long n,float *sx,long incx,float *sy,long incy)
|
|
{
|
|
long i,ix,iy,m,mp1;
|
|
float sdot,stemp;
|
|
stemp = sdot = 0.0;
|
|
if(n <= 0) return sdot;
|
|
if(incx == 1 && incy == 1) goto S20;
|
|
ix = iy = 1;
|
|
if(incx < 0) ix = (-n+1)*incx+1;
|
|
if(incy < 0) iy = (-n+1)*incy+1;
|
|
for(i=1; i<=n; i++) {
|
|
stemp += (*(sx+ix-1)**(sy+iy-1));
|
|
ix += incx;
|
|
iy += incy;
|
|
}
|
|
sdot = stemp;
|
|
return sdot;
|
|
S20:
|
|
m = n % 5L;
|
|
if(m == 0) goto S40;
|
|
for(i=0; i<m; i++) stemp += (*(sx+i)**(sy+i));
|
|
if(n < 5) goto S60;
|
|
S40:
|
|
mp1 = m+1;
|
|
for(i=mp1; i<=n; i+=5) stemp += (*(sx+i-1)**(sy+i-1)+*(sx+i)**(sy+i)+*(sx+i
|
|
+1)**(sy+i+1)+*(sx+i+2)**(sy+i+2)+*(sx+i+3)**(sy+i+3));
|
|
S60:
|
|
sdot = stemp;
|
|
return sdot;
|
|
}
|
|
void spofa(float *a,long lda,long n,long *info)
|
|
/*
|
|
SPOFA FACTORS A REAL SYMMETRIC POSITIVE DEFINITE MATRIX.
|
|
SPOFA IS USUALLY CALLED BY SPOCO, BUT IT CAN BE CALLED
|
|
DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
|
|
(TIME FOR SPOCO) = (1 + 18/N)*(TIME FOR SPOFA) .
|
|
ON ENTRY
|
|
A REAL(LDA, N)
|
|
THE SYMMETRIC MATRIX TO BE FACTORED. ONLY THE
|
|
DIAGONAL AND UPPER TRIANGLE ARE USED.
|
|
LDA INTEGER
|
|
THE LEADING DIMENSION OF THE ARRAY A .
|
|
N INTEGER
|
|
THE ORDER OF THE MATRIX A .
|
|
ON RETURN
|
|
A AN UPPER TRIANGULAR MATRIX R SO THAT A = TRANS(R)*R
|
|
WHERE TRANS(R) IS THE TRANSPOSE.
|
|
THE STRICT LOWER TRIANGLE IS UNALTERED.
|
|
IF INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
|
|
INFO INTEGER
|
|
= 0 FOR NORMAL RETURN.
|
|
= K SIGNALS AN ERROR CONDITION. THE LEADING MINOR
|
|
OF ORDER K IS NOT POSITIVE DEFINITE.
|
|
LINPACK. THIS VERSION DATED 08/14/78 .
|
|
CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
|
|
SUBROUTINES AND FUNCTIONS
|
|
BLAS SDOT
|
|
FORTRAN SQRT
|
|
INTERNAL VARIABLES
|
|
*/
|
|
{
|
|
extern float sdot(long n,float *sx,long incx,float *sy,long incy);
|
|
long j,jm1,k;
|
|
float t,s;
|
|
/*
|
|
BEGIN BLOCK WITH ...EXITS TO 40
|
|
*/
|
|
for(j=1; j<=n; j++) {
|
|
*info = j;
|
|
s = 0.0;
|
|
jm1 = j-1;
|
|
if(jm1 < 1) goto S20;
|
|
for(k=0; k<jm1; k++) {
|
|
t = *(a+k+(j-1)*lda)-sdot(k,(a+k*lda),1L,(a+(j-1)*lda),1L);
|
|
t /= *(a+k+k*lda);
|
|
*(a+k+(j-1)*lda) = t;
|
|
s += (t*t);
|
|
}
|
|
S20:
|
|
s = *(a+j-1+(j-1)*lda)-s;
|
|
/*
|
|
......EXIT
|
|
*/
|
|
if(s <= 0.0) goto S40;
|
|
*(a+j-1+(j-1)*lda) = sqrt(s);
|
|
}
|
|
*info = 0;
|
|
S40:
|
|
return;
|
|
}
|