Input matrix is assumed to have more rows than columns */
| 252 | |
| 253 | /* Input matrix is assumed to have more rows than columns */ |
| 254 | void pinv(imatrix& Ai,const imatrix& B) |
| 255 | { |
| 256 | imatrix A=B; |
| 257 | int nr=A.size(),nc=A[0].size(); |
| 258 | ident(Ai,nr); |
| 259 | |
| 260 | int r=0,c=0; |
| 261 | bool flag=true; |
| 262 | |
| 263 | bool verbose = OnlineOptions::singleton.has_option("verbose_he_setup"); |
| 264 | |
| 265 | if (verbose) |
| 266 | cout << "In inverse " << nr << " x " << nc << endl; |
| 267 | while (flag) |
| 268 | { //cout << "Inv step r=" << r << " c=" << c << endl; |
| 269 | if ((c%100)==0 and verbose) |
| 270 | cout << "Inv: " << c << " out of " << nc << endl; |
| 271 | // Find pivot |
| 272 | int k=r; |
| 273 | while (A[k][c]==0) { k++; } |
| 274 | // Swap rows if needed |
| 275 | if (k!=r) |
| 276 | { |
| 277 | A[r].swap(A[k]); |
| 278 | Ai[r].swap(Ai[k]); |
| 279 | } |
| 280 | // Kill off all rows above and below with a one in this position |
| 281 | for (k=0; k<nr; k++) |
| 282 | { if (k!=r && A[k][c]==1) |
| 283 | { // Only have to go from c onwards as rest are done |
| 284 | A[k].add(A[r]); |
| 285 | // Need to do all cols here |
| 286 | Ai[k].add(Ai[r]); |
| 287 | } |
| 288 | } |
| 289 | r++; c++; |
| 290 | if (c==nc) { flag=false; } |
| 291 | } |
| 292 | } |
| 293 | |
| 294 | |
| 295 | bool imatrix::operator!=(const imatrix& other) const |