PLearn 0.1
SparseMatrix.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PLearn (A C++ Machine Learning Library)
00004 // Copyright (C) 1999-2002 Yoshua Bengio
00005 //
00006 
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 // 
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 // 
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 // 
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 // 
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 
00036 #include "SparseMatrix.h"
00037 
00038 namespace PLearn {
00039 using namespace std;
00040 
00041 void SparseMatrix::resize(int nbrows,int n_columns,int n_non_zero)
00042 {
00043     n_rows=nbrows;
00044     beginRow.resize(n_columns);
00045     endRow.resize(n_columns);
00046     row.resize(n_non_zero);
00047     values.resize(n_non_zero);
00048 }
00049 
00050 // load SparseMatrix from file in ascii Harwell-Boeing Fortran format:
00051 // 4-line header, followed by beginRow, row, and values.
00052 void SparseMatrix::loadFortran(const char* filename)
00053 {
00054     FILE* fp=fopen(filename,"r");
00055     if (!fp)
00056         PLERROR("SparseMatrix::loadFortran, can't open file %s\n",filename);
00057     int n_cols,n_nonzero;
00058     fscanf(fp,"%*72c%*s%*s%*s%d%d%d%*d",&n_rows, &n_cols, &n_nonzero); 
00059     fscanf(fp,"%*s %*s %*s %*s"); // skip some format infos
00060     beginRow.resize(n_cols);
00061     endRow.resize(n_cols);
00062     values.resize(n_nonzero);
00063     row.resize(n_nonzero);
00064     real *brow = beginRow.data();
00065     real *erow = endRow.data();
00066     real *v = values.data();
00067     real *r = row.data();
00068     int i;
00069     for (i = 0; i < n_cols; i++) 
00070     { 
00071 #ifdef USEFLOAT
00072         fscanf(fp, "%f", &brow[i]); 
00073 #endif
00074 #ifdef USEDOUBLE
00075         fscanf(fp, "%lf", &brow[i]); 
00076 #endif
00077         brow[i]--; 
00078         if (i>0) erow[i-1]=brow[i]-1;
00079     }
00080     erow[n_cols-1]=n_nonzero-1;
00081     fscanf(fp,"%d",&i);
00082     if (i-1!=n_nonzero) 
00083         PLERROR("SparseMatrix::loadFortran, inconsistent nnz %d vs %d",
00084                 n_nonzero,i);
00085     for (i=0;i<n_nonzero;i++)
00086     {
00087 #ifdef USEFLOAT
00088         fscanf(fp, "%f", &r[i]); 
00089 #endif
00090 #ifdef USEDOUBLE
00091         fscanf(fp, "%lf", &r[i]); 
00092 #endif
00093         r[i]--;
00094     }
00095     for (i=0;i<n_nonzero;i++)
00096 #ifdef USEFLOAT
00097         fscanf(fp, "%f", &v[i]); 
00098 #endif
00099 #ifdef USEDOUBLE
00100     fscanf(fp, "%lf", &v[i]); 
00101 #endif
00102 }
00103 // save SparseMatrix from file in ascii Harwell-Boeing Fortran format:
00104 // 4-line header, followed by beginRow, row, and values.
00105 void SparseMatrix::saveFortran(const char* filename)
00106 {
00107     FILE* fp=fopen(filename,"w");
00108     if (!fp)
00109         PLERROR("SparseMatrix::saveFortran, can't open file %s\n",filename);
00110     int n_nonzero=values.length(), n_cols = endRow.length();
00111     fprintf(fp,"%72s%8s\n#\nrra %d %d %d 0\n","SparseMatrix         ",
00112             filename,
00113             n_rows, n_cols , n_nonzero);
00114     fprintf(fp,"          (10i8)          (10i8)            (8f10.3)            (8f10.3)\n");
00115     real *brow = beginRow.data();
00116     real *v = values.data();
00117     real *r = row.data();
00118     int i;
00119     for (i = 0; i < n_cols; i++) 
00120         //fprintf(fp, "%8d", (int)(brow[i]+1)); 
00121         fprintf(fp, "%7d ", (int)(brow[i]+1)); 
00122     fprintf(fp,"%8d\n",values.length()+1);
00123     for (i=0;i<n_nonzero;i++)
00124         fprintf(fp,"%7d ",(int)(r[i]+1));
00125     fprintf(fp,"\n");
00126     for (i=0;i<n_nonzero;i++)
00127         fprintf(fp,"%9f ",v[i]);
00128     fprintf(fp,"\n");
00129     fclose(fp);
00130 }
00131 
00132 // convert to the equivalent full matrix
00133 Mat SparseMatrix::toMat()
00134 {
00135     int n_cols = beginRow.length();
00136     Mat A(n_rows,n_cols);
00137     real* r=row.data();
00138     real* v=values.data();
00139     for (int c=0;c<n_cols;c++)
00140     {
00141         real* Ac = &A(0,c);
00142         int e = (int)endRow[c];
00143         for (int k=(int)beginRow[c];k<=e;k++)
00144             Ac[n_cols*(int)r[k]]=v[k];
00145     }
00146     return A;
00147 }
00148 
00149 SparseMatrix::SparseMatrix(Mat A)
00150     : n_rows(A.length()), beginRow(A.width()), endRow(A.width())
00151 {
00152     int n_nonzero=0;
00153     for (int i=0;i<n_rows;i++)
00154     {
00155         real* Ai=A[i];
00156         for (int j=0;j<A.width();j++)
00157             if (Ai[j]!=0) n_nonzero++;
00158     }
00159     row.resize(n_nonzero);
00160     values.resize(n_nonzero);
00161     int mod = A.mod();
00162     int k=0;
00163     real* v=values.data();
00164     real* r=row.data();
00165     real* b=beginRow.data();
00166     real* e=endRow.data();
00167     for (int j=0;j<A.width();j++)
00168     {
00169         real* Aij = &A(0,j);
00170         b[j] = k;
00171         for (int i=0;i<n_rows;i++,Aij+=mod)
00172             if (*Aij!=0)
00173             {
00174                 r[k] = i;
00175                 v[k] = *Aij;
00176                 k++;
00177             }
00178         e[j] = k-1;
00179     }
00180 }
00181 
00182 void SparseMatrix::product(const Vec& x, Vec& y)
00183 {
00184     // y[i] = sum_j A[i,j] x[j]
00185     if (y.length()!=n_rows || x.length()!=beginRow.length())
00186         PLERROR("SparseMatrix(%d,%d)::product(x(%d) -> y(%d)): dimensions don't match",
00187                 n_rows,beginRow.length(),x.length(),y.length());
00188     y.clear();
00189     real* y_=y.data();
00190     real* x_=x.data();
00191     real* A_=values.data();
00192     // loop over columns of A, accumulating in y
00193     for (int j=0;j<beginRow.length();j++)
00194     {
00195         real xj = x_[j];
00196         for (int k=(int)beginRow[j];k<=endRow[j];k++)
00197         {
00198             int i=(int)row[k];
00199             y_[i] += A_[k] * xj;
00200         }
00201     }
00202 }
00203 
00204 void SparseMatrix::diag(Vec& d)
00205 {
00206     real* d_ = d.data();
00207     real* A_ = values.data();
00208     int k;
00209     for (int j=0;j<beginRow.length();j++)
00210     {
00211         int end=int(endRow[j]);
00212         for (k=(int)beginRow[j];k<=end && int(row[k])!=j;k++);
00213         if (k<=end)
00214             d_[j]=A_[k];
00215         else
00216             d_[j]=0;
00217     }
00218 }
00219 
00220 void SparseMatrix::diagonalOfSquare(Vec& d)
00221 {
00222     real* d_ = d.data();
00223     real* A_ = values.data();
00224     int k;
00225     for (int j=0;j<beginRow.length();j++)
00226     {
00227         int end=int(endRow[j]);
00228         real sum2=0;
00229         for (k=(int)beginRow[j];k<=end;k++)
00230             sum2 += A_[k]*A_[k];
00231         d_[j]=sum2;
00232     }
00233 }
00234 
00235 // return dot product of i-th row with vector v
00236 real SparseMatrix::dotRow(int i, Vec v)
00237 {
00238     PLERROR("SparseMatrix is not appropriate to perform dotRow operations");
00239     return 0;
00240 }
00241 
00242 // return dot product of j-th column with vector v
00243 real SparseMatrix::dotColumn(int j, Vec v)
00244 {
00245 #ifdef BOUNDCHECK
00246     if (v.length()!=length())
00247         PLERROR("SparseMatrix::dotColumn(%d,v), v.length_=%d != matrix length=%d",
00248                 j,v.length(),length());
00249 #endif
00250     real s=0;
00251     real* v_=v.data();
00252     real* A_=values.data();
00253     for (int k=int(beginRow[j]);k<=int(endRow[j]);k++)
00254         s += A_[k] * v_[int(row[k])];
00255     return s;
00256 }
00257 
00258 // add two sparse matrices (of same dimensions)
00259 SparseMatrix operator+(const SparseMatrix& A, const SparseMatrix& B)
00260 {
00261     int n_rows = A.n_rows;
00262     int n_columns = A.beginRow.length();
00263     if (n_rows != B.n_rows)
00264         PLERROR("SparseMatrix(%d,%d)+SparseMatrix(%d,%d): both should have same dimensions",
00265                 n_rows,A.beginRow.length(),B.n_rows,B.beginRow.length());
00266     int n_non_zero = A.row.length()+B.row.length(); // THIS IS AN UPPER BOUND ON ACTUAL n_non_zero
00267     SparseMatrix C(n_rows,n_columns,n_non_zero);
00268 
00269     int n_actual_non_zero=0;
00270     // the data is stored column-wise
00271     Vec column(n_rows);
00272     real* v=column.data();
00273     for (int j=0;j<n_columns;j++)
00274     {
00275         column.clear();
00276         for (int i=(int)A.beginRow[j];i<=A.endRow[j];i++)
00277             v[(int)A.row[i]]=A.values[i];
00278         for (int i=(int)B.beginRow[j];i<=B.endRow[j];i++)
00279             v[(int)B.row[i]]+=B.values[i];
00280         C.beginRow[j]=n_actual_non_zero;
00281         for (int i=0;i<n_rows;i++)
00282             if (v[i]!=0)
00283             {
00284                 C.row[n_actual_non_zero]=i;
00285                 C.values[n_actual_non_zero]=v[i];
00286                 n_actual_non_zero++;
00287             }
00288         C.endRow[j]=n_actual_non_zero-1;
00289 
00290     }
00291     C.row.resize(n_actual_non_zero);
00292     C.values.resize(n_actual_non_zero);
00293     return C;
00294 }
00295 
00296 // add a bunch of sparse matrices and return result
00297 SparseMatrix add(Array<SparseMatrix>& matrices)
00298 {
00299     int n_mat = matrices.size();
00300     if (n_mat<1) PLERROR("add(Array<SparseMatrix>) argument is empty!");
00301     int n_rows = matrices[0].n_rows;
00302     int n_columns = matrices[0].beginRow.length();
00303     for (int i=1;i<n_mat;i++)
00304         if (n_rows != matrices[i].n_rows)
00305             PLERROR("add(SparseMatrix(%d,%d)+SparseMatrix(%d,%d): both should have same dimensions",
00306                     n_rows,n_columns,matrices[i].n_rows,matrices[i].beginRow.length());
00307     int n_non_zero = 0;
00308     for (int i=0;i<n_mat;i++)
00309         n_non_zero+=matrices[i].row.length(); // UPPER BOUND ON ACTUAL n_non_zero
00310 
00311     SparseMatrix C(n_rows,n_columns,n_non_zero);
00312 
00313     int n_actual_non_zero=0;
00314     // the data is stored column-wise
00315     Vec column(n_rows);
00316     real* v=column.data();
00317     for (int j=0;j<n_columns;j++)
00318     {
00319         column.clear();
00320         for (int k=0;k<n_mat;k++)
00321             for (int i=(int)matrices[k].beginRow[j];i<=(int)matrices[k].endRow[j];i++)
00322                 v[(int)matrices[k].row[i]]+=matrices[k].values[i];
00323         C.beginRow[j]=n_actual_non_zero;
00324         for (int i=0;i<n_rows;i++)
00325             if (v[i]!=0)
00326             {
00327                 C.row[n_actual_non_zero]=i;
00328                 C.values[n_actual_non_zero]=v[i];
00329                 n_actual_non_zero++;
00330             }
00331         C.endRow[j]=n_actual_non_zero-1;
00332 
00333     }
00334     C.row.resize(n_actual_non_zero);
00335     C.values.resize(n_actual_non_zero);
00336     return C;
00337 }
00338 
00339 } // end of namespace PLearn
00340 
00341 
00342 /*
00343   Local Variables:
00344   mode:c++
00345   c-basic-offset:4
00346   c-file-style:"stroustrup"
00347   c-file-offsets:((innamespace . 0)(inline-open . 0))
00348   indent-tabs-mode:nil
00349   fill-column:79
00350   End:
00351 */
00352 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines