sparse.cpp 3.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177
  1. /******* functions, needed for sparse matrices *******/
  2. #include <math.h>
  3. #include <iostream.h>
  4. #include <stdlib.h>
  5. #include <iomanip.h>
  6. #include "complex.h"
  7. #include "vector.h"
  8. #include "matrix.h"
  9. #include "cmplxvec.h"
  10. #include "cmplxmat.h"
  11. #include "ivectorl.h"
  12. #include "sparse.h"
  13. //#include "stdafx.h"
  14. void sprsin_d( Matrix& a, double thresh, Vector& sa, IVectorl& ija ) {
  15. const int ns = a.dim_i();
  16. int i, j, k, nh;
  17. nh = ns + 1;
  18. for( i = 0; i < ns; i++ ) {
  19. for( j = 0; j < ns; j++ ) {
  20. if( fabs( a(i,j) ) > thresh && i != j )
  21. nh += 1;
  22. }
  23. }
  24. sa.resize( nh );
  25. ija.resize( nh );
  26. for( j = 0; j < ns; j++ )
  27. sa[j] = a(j,j);
  28. ija[0] = ns + 1;
  29. k = ns;
  30. for( i = 0; i < ns; i++ ) {
  31. for( j = 0; j < ns; j++ ) {
  32. if( fabs( a(i,j) ) > thresh && i != j ) {
  33. k += 1;
  34. sa[k] = a(i,j);
  35. ija[k] = j;
  36. }
  37. }
  38. ija[i+1] = k + 1;
  39. }
  40. a.resize(0,0);
  41. }
  42. void sprsax_d( Vector& sa, IVectorl& ija, Vector& x, Vector& b ) {
  43. const int ns = x.dim();
  44. int i, k;
  45. b.resize(ns);
  46. if( ija[0] != ns + 1 ) {
  47. cerr << " Something is wrong in sprsax_d !!! " << endl;
  48. exit(1);
  49. }
  50. for( i = 0; i < ns; i++ ) {
  51. b[i] = sa[i] * x[i]; //cerr << i << endl;
  52. for( k = ija[i]; k <= ija[i+1]-1; k++ ) {
  53. //cerr << k << endl;
  54. b[i] += sa[k] * x[ija[k]];
  55. }
  56. }
  57. }
  58. void sprsin_c( CmplxMatrix& a, double thresh, CmplxVector& sa, IVectorl& ija ) {
  59. const int ns = a.dim_i();
  60. int i, j, k, nh;
  61. nh = ns + 1;
  62. for( i = 0; i < ns; i++ ) {
  63. for( j = 0; j < ns; j++ ) {
  64. if( cabs( a(i,j) ) > thresh && i != j )
  65. nh += 1;
  66. }
  67. }
  68. sa.resize( nh );
  69. ija.resize( nh );
  70. for( j = 0; j < ns; j++ )
  71. sa[j] = a(j,j);
  72. ija[0] = ns + 1;
  73. k = ns;
  74. for( i = 0; i < ns; i++ ) {
  75. for( j = 0; j < ns; j++ ) {
  76. if( cabs( a(i,j) ) > thresh && i != j ) {
  77. k += 1;
  78. sa[k] = a(i,j);
  79. ija[k] = j;
  80. }
  81. }
  82. ija[i+1] = k + 1;
  83. }
  84. a.resize(0,0);
  85. }
  86. void sprsax_c( CmplxVector& sa, IVectorl& ija, CmplxVector& x, CmplxVector& b ) {
  87. const int ns = x.dim();
  88. int i, k;
  89. b.resize(ns);
  90. if( ija[0] != ns + 1 ) {
  91. cerr << " Something is wrong in sprsax_c !!! " << endl;
  92. exit(1);
  93. }
  94. for( i = 0; i < ns; i++ ) {
  95. b[i] = sa[i] * x[i];
  96. for( k = ija[i]; k <= ija[i+1]-1; k++ ) {
  97. b[i] += sa[k] * x[ija[k]];
  98. }
  99. }
  100. }
  101. Vector sprsax_d_v( Vector& sa, IVectorl& ija, Vector& x ) {
  102. Vector res;
  103. sprsax_d( sa, ija, x, res );
  104. return( res );
  105. }
  106. CmplxVector sprsax_c_v( CmplxVector& sa, IVectorl& ija, CmplxVector& x ) {
  107. CmplxVector res;
  108. sprsax_c( sa, ija, x, res );
  109. return( res );
  110. }