00001 #include "csv.h"
00002 #include <stdlib.h>
00003 #include <stdio.h>
00004 #include <math.h>
00005
00006 compvect::compvect ()
00007 {
00008 n=0; nz=0;
00009 a=NULL;
00010 ind=NULL;
00011
00012
00013 limit=1.0e-15;
00014
00015 ordering=0;
00016 }
00017
00018 compvect::compvect (long *ii,long nonz)
00019 {
00020 long i;
00021
00022 n=0; nz=nonz;
00023 a=new double [nz];
00024 ind=new long [nz];
00025
00026 for (i=0;i<nz;i++){
00027 ind[i]=ii[i];
00028 a[i]=0.0;
00029 }
00030
00031
00032 limit=1.0e-15;
00033
00034 ordering=0;
00035 }
00036
00037 compvect::~compvect ()
00038 {
00039 delete [] a;
00040 delete [] ind;
00041 }
00042
00043
00044
00045
00046
00047
00048
00049
00050 void compvect::copy (compvect *bc)
00051 {
00052 long i;
00053
00054 bc->~compvect ();
00055
00056 bc->n = n;
00057 bc->nz = nz;
00058 bc->ind = new long [nz];
00059 bc->a = new double [nz];
00060
00061 for (i=0;i<nz;i++){
00062 bc->ind[i] = ind[i];
00063 bc->a[i] = a[i];
00064 }
00065
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 void compvect::setup_vector (double *b,long nc)
00077 {
00078 long i;
00079
00080
00081 n=nc;
00082
00083
00084 nz=0;
00085 for (i=0;i<n;i++){
00086 if (fabs(b[i])>limit)
00087 nz++;
00088 }
00089
00090 if (a!=NULL){
00091 delete [] a;
00092 }
00093 a = new double [nz];
00094 if (ind!=NULL){
00095 delete [] ind;
00096 }
00097 ind = new long [nz];
00098
00099 nz=0;
00100 for (i=0;i<n;i++){
00101 if (fabs(b[i])>limit){
00102 ind[nz]=i;
00103 a[nz]=b[i];
00104 nz++;
00105 }
00106 }
00107
00108
00109 ordering=1;
00110 }
00111
00112
00113
00114
00115
00116
00117 void compvect::test_ordering ()
00118 {
00119 long i;
00120
00121 for (i=1;i<nz;i++){
00122 if (ind[i]<=ind[i-1])
00123 break;
00124 }
00125
00126 ordering=0;
00127 }
00128
00129
00130
00131
00132
00133
00134
00135 void compvect::reorder_components ()
00136 {
00137 long i,j,k,l,min;
00138 double s;
00139
00140 for (i=0;i<nz;i++){
00141 min=ind[i];
00142 for (j=i;j<nz;j++){
00143 if (min>ind[j]){
00144 min=ind[j];
00145 k=j;
00146 }
00147 }
00148 l=ind[i];
00149 ind[i]=k;
00150 ind[k]=l;
00151
00152 s=a[i];
00153 a[i]=min;
00154 a[k]=s;
00155 }
00156
00157
00158 ordering=1;
00159 }
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 void compvect::axpy (compvect *x,double aa)
00170 {
00171 long i,j,k;
00172 compvect bc;
00173
00174
00175 copy (&bc);
00176
00177 i=0; j=0; k=0;
00178 do{
00179 if (ind[i] == x->ind[j]){
00180 i++; j++; k++;
00181 }
00182 else{
00183 if (ind[i] < x->ind[j]){
00184 i++;
00185 }
00186 else{
00187 j++;
00188 }
00189 }
00190 }while(i<nz && j<x->nz);
00191
00192 nz=k;
00193 delete [] ind;
00194 delete [] a;
00195 ind = new long [nz];
00196 a = new double [nz];
00197
00198 i=0; j=0; k=0;
00199 do{
00200 if (bc.ind[i] == x->ind[j]){
00201 a[k]=bc.a[i]+x->a[j]*aa;
00202 ind[k]=bc.ind[i];
00203 i++; j++; k++;
00204 }
00205 else{
00206 if (bc.ind[i] < x->ind[j]){
00207 i++;
00208 }
00209 else{
00210 j++;
00211 }
00212 }
00213 }while(i<bc.nz && j<x->nz);
00214
00215 }
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 void compvect::axpy_known (compvect *x,double aa)
00226 {
00227 long i;
00228
00229 for (i=0;i<nz;i++){
00230 a[i] += aa*x->a[i];
00231 }
00232 }
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 double compvect::dotprod (compvect *cv)
00244 {
00245 long j,k;
00246 long min1,max1,min2,max2,nc;
00247 long *ind2;
00248 double s,*a2;
00249
00250 if (cv->ordering == 0){
00251 fprintf (stderr,"\n\n compressed sparse vector is not ordered increasingly in function comprow::ss (file %s, line %d)\n",__FILE__,__LINE__);
00252 }
00253 if (ordering == 0){
00254 fprintf (stderr,"\n\n compressed sparse vector is not ordered increasingly in function comprow::ss (file %s, line %d)\n",__FILE__,__LINE__);
00255 }
00256
00257 ind2=cv->ind;
00258 a2=cv->a;
00259
00260
00261 nc=cv->nz;
00262
00263 min2=ind2[0];
00264
00265 max2=ind2[nc-1];
00266
00267
00268 min1=ind[0];
00269
00270 max1=ind[nz-1];
00271
00272 if (max1<min2)
00273 return s=0.0;
00274 if (max2<min1)
00275 return s=0.0;
00276
00277 j=0; k=0;
00278 s=0.0;
00279 do{
00280 if (ind[j]==ind2[k]){
00281 s+=a[j]*a2[k];
00282 j++; k++;
00283 }
00284 else{
00285 if (ind[j]<ind[k]){
00286 j++;
00287 }
00288 else{
00289 k++;
00290 }
00291 }
00292 }while (j<nz && k<nc);
00293
00294 return s;
00295 }