/[escript]/branches/lapack2681/escript/src/DataExpanded.cpp
ViewVC logotype

Diff of /branches/lapack2681/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2739 by jfenwick, Mon Nov 9 06:25:31 2009 UTC revision 2740 by jfenwick, Tue Nov 10 06:48:24 2009 UTC
# Line 567  DataExpanded::matrixInverse(DataAbstract Line 567  DataExpanded::matrixInverse(DataAbstract
567      throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
568    }    }
569    
570      if (getRank()!=2)
571      {
572        throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
573    
574      }
575    int  sampleNo;    int  sampleNo;
576    const int numdpps=getNumDPPSample();    const int numdpps=getNumDPPSample();
577    const int numSamples = getNumSamples();    const int numSamples = getNumSamples();
578    const ValueType& vec=m_data.getData();    const ValueType& vec=m_data.getData();
579    #pragma omp parallel for private(sampleNo,p) schedule(static)    int errorcode=0;
580    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    #pragma omp parallel private(sampleNo) reduction(MAX:errorcode)
581              // not sure I like all those virtual calls    {
582      DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);       int* p=new int[getShape()[0]];
583      DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps);       scoped_ptr<int> piv(p);
584         #pragma omp parallel for private(sampleNo) schedule(static)
585         for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
586         {
587                // not sure I like all those virtual calls to getPointOffset
588            DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
589            int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, p);
590        errorcode=(res>errorcode)?res:errorcode;
591         }
592      }
593      if (errorcode)
594      {
595        DataMaths::matrixInverseError(errorcode);   // throws exceptions
596    }    }
597  }  }
598    

Legend:
Removed from v.2739  
changed lines
  Added in v.2740

  ViewVC Help
Powered by ViewVC 1.1.26