ireal acd_rdom_inner(std::vector<RDOM *> updt, std::vector <RDOM*> old, bool fwd){
    try{
      if (updt.size() != old.size()){
	RVLException e;
	e<<"Error:acd_rdom_inner -- two input vector have different sizes.\n";
	throw e;
      }
      ireal dot=REAL_ZERO;
      int n = updt.size();
        if (fwd==true) {
            ireal axy, ip;
            if (n == 2) {
                RARR *upd0_n = &(updt[1]->_s[D_UC]);
                RARR *ucd0_n = &(updt[1]->_s[D_UP]);
                RARR *upb_o = &(old[1]->_s[D_UP]);
                RARR *ucb_o = &(old[1]->_s[D_UC]);
                
                ra_a_inner(upd0_n, upb_o, &axy);
                ra_a_inner(ucd0_n, ucb_o, &ip);
                
                axy = axy + ip;
            }
            if (n == 4) {
                RARR *upd0_n = &(updt[2]->_s[D_UC]);
                RARR *ucd0_n = &(updt[2]->_s[D_UP]);
                RARR *upb_o = &(old[2]->_s[D_UP]);
                RARR *ucb_o = &(old[2]->_s[D_UC]);
                
                ra_a_inner(upd0_n, upb_o, &axy);
                ra_a_inner(ucd0_n, ucb_o, &ip);
                
                axy = axy + ip;
                
                RARR *updd_n = &(updt[n-1]->_s[D_UC]);
                RARR *ucdd_n = &(updt[n-1]->_s[D_UP]);
                RARR *updb_o = &(old[n-1]->_s[D_UP]);
                RARR *ucdb_o = &(old[n-1]->_s[D_UC]);
                
                ireal ip1;
                ra_a_inner(updd_n, updb_o, &ip1);
                ra_a_inner(ucdd_n, ucdb_o, &ip);
                
                axy = axy + ip + ip1;
            }
            dot=axy;
        }
        else if (fwd == false){
            ireal xaty, ip, ip1;
            if (n==2) {
                
                RARR *ucd0_o  = &(old[1]->_s[D_UC]);
                RARR *upd0_o  = &(old[1]->_s[D_UP]);
                RARR *csqd0_o = &(old[1]->_s[D_CSQ]);
                RARR *ucb_n   = &(updt[1]->_s[D_UC]);
                RARR *upb_n   = &(updt[1]->_s[D_UP]);
                RARR *csqb_n  = &(updt[1]->_s[D_CSQ]);
                
                ra_a_inner(upd0_o, upb_n, &xaty);
                ra_a_inner(csqd0_o, csqb_n, &ip);
                ra_a_inner(ucd0_o, ucb_n, &ip1);
                xaty = xaty + ip + ip1;
            }
            if (n == 4) {
                RARR *ucd0_o  = &(old[2]->_s[D_UC]);
                RARR *upd0_o  = &(old[2]->_s[D_UP]);
                RARR *csqd0_o = &(old[2]->_s[D_CSQ]);
                RARR *ucb_n   = &(updt[2]->_s[D_UC]);
                RARR *upb_n   = &(updt[2]->_s[D_UP]);
                RARR *csqb_n  = &(updt[2]->_s[D_CSQ]);
                
                ra_a_inner(upd0_o, upb_n, &xaty);
                ra_a_inner(csqd0_o, csqb_n, &ip);
                ra_a_inner(ucd0_o, ucb_n, &ip1);
                xaty = xaty + ip + ip1;
                
                RARR *updd_o  = &(old[n-1]->_s[D_UP]);
                RARR *ucdd_o  = &(old[n-1]->_s[D_UC]);
                RARR *ucdb_n  = &(updt[n-1]->_s[D_UC]);
                RARR *updb_n  = &(updt[n-1]->_s[D_UP]);
                
                ra_a_inner(ucdd_o, ucdb_n, &ip);
                ra_a_inner(updd_o, updb_n, &ip1);
                xaty = xaty + ip + ip1;
            }
            dot=xaty;
        }

/*      if (fwd==true) {
	RARR *upd_n = &(updt[n-1]->_s[D_UC]);
	RARR *ucd_n = &(updt[n-1]->_s[D_UP]);
	RARR *upb_o = &(old[n-1]->_s[D_UP]);
	RARR *ucb_o = &(old[n-1]->_s[D_UC]);
	ireal axy, ip;
	ra_a_inner(upd_n, upb_o, &axy);
	ra_a_inner(ucd_n, ucb_o, &ip);
	axy = axy + ip;
            
	return axy;
      }
      else if (fwd == false){
	RARR *updd_o  = &(old[n-1]->_s[D_UP]);
	RARR *ucdd_o  = &(old[n-1]->_s[D_UC]);
	RARR *csqdd_o = &(old[n-1]->_s[D_CSQ]);
	RARR *ucdb_n  = &(updt[n-1]->_s[D_UC]);
	RARR *updb_n  = &(updt[n-1]->_s[D_UP]);
	RARR *csqdb_n = &(updt[n-1]->_s[D_CSQ]);
	ireal xaty, ip, ip1;
	ra_a_inner(updd_o, updb_n, &xaty);
	ra_a_inner(csqdd_o, csqdb_n, &ip);
	ra_a_inner(ucdd_o, ucdb_n, &ip1);
	xaty = xaty + ip + ip1;
	if (n == 4) {
	  RARR *ucd0_o  = &(old[1]->_s[D_UC]);
	  RARR *csqd0_o = &(old[1]->_s[D_CSQ]);
	  RARR *ucb_n   = &(updt[1]->_s[D_UC]);
	  RARR *csqb_n  = &(updt[1]->_s[D_CSQ]);
	  ra_a_inner(ucd0_o, ucb_n, &ip);
	  ra_a_inner(csqd0_o, csqb_n, &ip1);
	  xaty = xaty + ip + ip1;
	}
	return xaty;
      }
*/
	return dot;
    }
    catch (RVLException & e) {
      e.write(cerr);
      exit(1);
    }
  }
    ireal acd_rdom_inner(std::vector<RDOM *> updt, std::vector <RDOM> old, bool fwd){
        try{
            if (updt.size() != old.size()){
                RVLException e;
                e<<"Error:acd_rdom_inner -- two input vector have different sizes.\n";
                throw e;
            }
	    ireal dot=REAL_ZERO;
            int n = updt.size();
            if (fwd==true) {
                ireal axy, ip;
                if (n == 2) {
                RARR *upd0_n = &(updt[1]->_s[D_UC]);
                RARR *ucd0_n = &(updt[1]->_s[D_UP]);
                RARR *upb_o = &(old[1]._s[D_UP]);
                RARR *ucb_o = &(old[1]._s[D_UC]);
                
                ra_a_inner(upd0_n, upb_o, &axy);
                ra_a_inner(ucd0_n, ucb_o, &ip);
                
                axy = axy + ip;
                }
                if (n == 4) {
                    RARR *upd0_n = &(updt[2]->_s[D_UC]);
                    RARR *ucd0_n = &(updt[2]->_s[D_UP]);
                    RARR *upb_o = &(old[2]._s[D_UP]);
                    RARR *ucb_o = &(old[2]._s[D_UC]);
                    
                    ra_a_inner(upd0_n, upb_o, &axy);
                    ra_a_inner(ucd0_n, ucb_o, &ip);
                    
                    axy = axy + ip;

                    RARR *updd_n = &(updt[n-1]->_s[D_UC]);
                    RARR *ucdd_n = &(updt[n-1]->_s[D_UP]);
                    RARR *updb_o = &(old[n-1]._s[D_UP]);
                    RARR *ucdb_o = &(old[n-1]._s[D_UC]);
                    
                    ireal ip1;
                    ra_a_inner(updd_n, updb_o, &ip1);
                    ra_a_inner(ucdd_n, ucdb_o, &ip);
                    
                    axy = axy + ip + ip1;
                }
                dot=axy;
            }
            else if (fwd == false){
                ireal xaty, ip, ip1;
                if (n==2) {
                
                RARR *ucd0_o  = &(old[1]._s[D_UC]);
                RARR *upd0_o  = &(old[1]._s[D_UP]);
                RARR *csqd0_o = &(old[1]._s[D_CSQ]);
                RARR *ucb_n   = &(updt[1]->_s[D_UC]);
                RARR *upb_n   = &(updt[1]->_s[D_UP]);
                RARR *csqb_n  = &(updt[1]->_s[D_CSQ]);
                
                ra_a_inner(upd0_o, upb_n, &xaty);
                ra_a_inner(csqd0_o, csqb_n, &ip);
                ra_a_inner(ucd0_o, ucb_n, &ip1);
                xaty = xaty + ip + ip1;
                }
                if (n == 4) {
                    RARR *ucd0_o  = &(old[2]._s[D_UC]);
                    RARR *upd0_o  = &(old[2]._s[D_UP]);
                    RARR *csqd0_o = &(old[2]._s[D_CSQ]);
                    RARR *ucb_n   = &(updt[2]->_s[D_UC]);
                    RARR *upb_n   = &(updt[2]->_s[D_UP]);
                    RARR *csqb_n  = &(updt[2]->_s[D_CSQ]);
                    
                    ra_a_inner(upd0_o, upb_n, &xaty);
                    ra_a_inner(csqd0_o, csqb_n, &ip);
                    ra_a_inner(ucd0_o, ucb_n, &ip1);
                    xaty = xaty + ip + ip1;
                    
                    RARR *updd_o  = &(old[n-1]._s[D_UP]);
                    RARR *ucdd_o  = &(old[n-1]._s[D_UC]);
                    RARR *ucdb_n  = &(updt[n-1]->_s[D_UC]);
                    RARR *updb_n  = &(updt[n-1]->_s[D_UP]);

                    ra_a_inner(ucdd_o, ucdb_n, &ip);
                    ra_a_inner(updd_o, updb_n, &ip1);
                    xaty = xaty + ip + ip1;
                }
                dot=xaty;
            }
	    return dot;
        }
        catch (RVLException & e) {
            e.write(cerr);
            exit(1);
        }
    }

  ireal acd_rdom_norm(std::vector<RDOM *> rd){
    int n=rd.size();
    RARR *upb = &(rd[n-1]->_s[D_UP]);
    RARR *ucb = &(rd[n-1]->_s[D_UC]);
        
    ireal norm, tmp;
    ra_a_inner(upb, upb, &norm);
    ra_a_inner(ucb, ucb, &tmp);
        
    norm = sqrt(norm + tmp);
    return norm;
  }
    
    ireal acd_rdom_inner_csq(RDOM *rd, int cflag){
        RARR *upb = &(rd->_s[D_UP]);
        RARR *ucb = &(rd->_s[D_UC]);

        ireal norm, tmp, tmp1;
        tmp1=0.0;
        ra_a_inner(upb, upb, &norm);
        ra_a_inner(ucb, ucb, &tmp);
        if (cflag==1) {
            RARR *csq = &(rd->_s[D_CSQ]);
            ra_a_inner(csq, csq, &tmp1);
        }
        norm = norm + tmp + tmp1;
        return norm;
    }
    ireal acd_rdom_inner_csq(RDOM rd, int cflag){
        RARR *upb = &(rd._s[D_UP]);
        RARR *ucb = &(rd._s[D_UC]);
        
        ireal norm, tmp, tmp1;
        tmp1=0.0;
        ra_a_inner(upb, upb, &norm);
        ra_a_inner(ucb, ucb, &tmp);
        if (cflag==1) {
            RARR *csq = &(rd._s[D_CSQ]);
            ra_a_inner(csq, csq, &tmp1);
        }
        norm = norm + tmp + tmp1;
        return norm;
    }
  ireal acd_rdom_norm(std::vector<RDOM> rd){
    try{
      int n=rd.size();
      RARR *upb = &(rd[n-1]._s[D_UP]);
      RARR *ucb = &(rd[n-1]._s[D_UC]);
            
      ireal norm, tmp;
      ra_a_inner(upb, upb, &norm);
      ra_a_inner(ucb, ucb, &tmp);
            
      norm = sqrt(norm + tmp);
      return norm;
    }
    catch (RVLException & e) {
      e.write(cerr);
      exit(1);
    }
  }