Mercurial > fem-fenics-eugenio

comparison src/assemble.cc @ 167:3895c8d4c066

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Improved assembling of the matrix.
author gedeone-octave <marcovass89@hotmail.it>
date Sun, 06 Oct 2013 22:23:40 +0100
parents 40ae9e5dfb93
children 67944f307560
comparison
equal deleted inserted replaced
166:40ae9e5dfb93 167:3895c8d4c066
95 } 95 }
96 else 96 else
97 error ("assemble: unknown argument type"); 97 error ("assemble: unknown argument type");
98 } 98 }
99 99
100 // It provides an estimation for excess of nnz elements 100 // It provides an upper boung for the nnz elements
101 boost::tuples::tuple<const std::size_t*, const std::size_t*, const double*, int> aa = A.data (); 101 boost::tuples::tuple<const std::size_t*, const std::size_t*, const double*, int> aa = A.data ();
102 int nnz = aa.get<3> (); 102 int nnz = aa.get<3> ();
103 std::size_t nr = A.size (0), nc = A.size (1); 103 std::size_t nr = A.size (0), nc = A.size (1);
104 std::vector<double> data_tmp; 104 std::vector<double> data_tmp;
105 std::vector<std::size_t> cidx_tmp; 105 std::vector<std::size_t> cidx_tmp;
106 106
107 octave_idx_type dims = nnz, nz = 0, ii = 0; 107 dim_vector dims (nnz, 1);
108 ColumnVector ridx (dims), cidx (dims), data (dims); 108 octave_idx_type nz = 0, ii = 0;
109 Array<octave_idx_type> ridx (dims, 0), cidx (dims, 0);
110 Array<double> data (dims, 0);
111
112 octave_idx_type* orow = ridx.fortran_vec ();
113 octave_idx_type* oc = cidx.fortran_vec ();
114 double* ov = data.fortran_vec ();
109 115
110 for (std::size_t i = 0; i < nr; ++i) 116 for (std::size_t i = 0; i < nr; ++i)
111 { 117 {
112 A.getrow (i, cidx_tmp, data_tmp); 118 A.getrow (i, cidx_tmp, data_tmp);
113 nz += cidx_tmp.size (); 119 nz += cidx_tmp.size ();
114 if (dims < nz)
115 {
116 dims = 1.2 * ((nr * nz) / (i + 1));;
117 ridx.resize (dims);
118 cidx.resize (dims);
119 data.resize (dims);
120 }
121 120
122 for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j) 121 for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j)
123 { 122 {
124 ridx.xelem (ii + j) = i + 1; 123 orow [ii + j] = i;
125 cidx.xelem (ii + j) = cidx_tmp [j] + 1; 124 oc [ii + j] = cidx_tmp [j];
126 data.xelem (ii + j) = data_tmp [j]; 125 ov [ii + j] = data_tmp [j];
127 } 126 }
128 127
129 ii = nz; 128 ii = nz;
130 } 129 }
131 130
132 ridx.resize (ii); 131 dims(0) = ii;
133 cidx.resize (ii); 132 ridx.resize (dims);
134 data.resize (ii); 133 cidx.resize (dims);
134 data.resize (dims);
135 135
136 SparseMatrix sm (data, ridx, cidx, nr, nc); 136 SparseMatrix sm (data, ridx, cidx, nr, nc);
137 retval(0) = sm; 137 retval(0) = sm;
138 } 138 }
139 139