psinv.c 1.2 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546
  1. /* psinv.c CCMATH mathematics library source code.
  2. *
  3. * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
  4. * This code may be redistributed under the terms of the GNU library
  5. * public license (LGPL). ( See the lgpl.license file for details.)
  6. * ------------------------------------------------------------------------
  7. */
  8. #include "ccmath.h"
  9. int psinv(double *v, int n)
  10. {
  11. double z, *p, *q, *r, *s, *t;
  12. int j, k;
  13. for (j = 0, p = v; j < n; ++j, p += n + 1) {
  14. for (q = v + j * n; q < p; ++q)
  15. *p -= *q * *q;
  16. if (*p <= 0.)
  17. return -1;
  18. *p = sqrt(*p);
  19. for (k = j + 1, q = p + n; k < n; ++k, q += n) {
  20. for (r = v + j * n, s = v + k * n, z = 0.; r < p;)
  21. z += *r++ * *s++;
  22. *q -= z;
  23. *q /= *p;
  24. }
  25. }
  26. trnm(v, n);
  27. for (j = 0, p = v; j < n; ++j, p += n + 1) {
  28. *p = 1. / *p;
  29. for (q = v + j, t = v; q < p; t += n + 1, q += n) {
  30. for (s = q, r = t, z = 0.; s < p; s += n)
  31. z -= *s * *r++;
  32. *q = z * *p;
  33. }
  34. }
  35. for (j = 0, p = v; j < n; ++j, p += n + 1) {
  36. for (q = v + j, t = p - j; q <= p; q += n) {
  37. for (k = j, r = p, s = q, z = 0.; k < n; ++k)
  38. z += *r++ * *s++;
  39. *t++ = (*q = z);
  40. }
  41. }
  42. return 0;
  43. }